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Abstract. This article reports on the confluence of two streams of research, 
one emanating from the fields of numerical analysis and scientific computa- 
tion, the other from topology and geometry. In it we consider the numerical 
discretization of partial differential equations that are related to differential 
complexes so that de Rham cohomology and Hodge theory are key tools for 
exploring the well-posedness of the continuous problem. The discretization 
methods we consider are finite element methods, in which a variational or 
weak formulation of the PDE problem is approximated by restricting the trial 
subspace to an appropriately constructed piecewise polynomial subspace. Af- 
ter a brief introduction to finite element methods, we develop an abstract 
Hilbert space framework for analyzing the stability and convergence of such 
discretizations. In this framework, the differential complex is represented by a 
complex of Hilbert spaces and stability is obtained by transferring Hodge the- 
oretic structures that ensure well-posedness of the continuous problem from 
the continuous level to the discrete. We show stable discretization is achieved 
if the finite element spaces satisfy two hypotheses: they can be arranged into a 
subcomplex of this Hilbert complex, and there exists a bounded cochain pro- 
jection from that complex to the subcomplex. In the next part of the paper, 
we consider the most canonical example of the abstract theory, in which the 
Hilbert complex is the de Rham complex of a domain in Euclidean space. We 
use the Koszul complex to construct two families of finite element differential 
forms, show that these can be arranged in subcomplexes of the de Rham com- 
plex in numerous ways, and for each construct a bounded cochain projection. 
The abstract theory therefore applies to give the stability and convergence 
of finite element approximations of the Hodge Laplacian. Other applications 
are considered as well, especially the elasticity complex and its application 
to the equations of elasticity. Background material is included to make the 
presentation self-contained for a variety of readers. 
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1. Introduction 



Numerical algorithms for the solution of partial differential equations are an es- 
sential tool of the modern world. They are applied in countless ways every day in 
problems as varied as the design of aircraft, prediction of climate, development of 
cardiac devices, and modeling of the financial system. Science, engineering, and 
technology depend not only on efficient and accurate algorithms for approximat- 
ing the solutions of a vast and diverse array of differential equations which arise 
in applications, but also on mathematical analysis of the behavior of computed 
solutions, in order to validate the computations, determine the ranges of apph- 
cability, compare different algorithms, and point the way to improved numerical 
methods. Given a partial differential equation (PDE) problem, a numerical algo- 
rithm approximates the solution by the solution of a finite dimensional problem 
which can be implemented and solved on a computer. This discretization depends 
on a parameter — representing, for example, a grid spacing, mesh size, or time step — 
which can be adjusted to obtain a more accurate approximate solution, at the cost 
of a larger finite dimensional problem. The mathematical analysis of the algorithm 
aims to describe the relationship between the true solution and the numerical so- 
lution as the discretization parameter is varied. For example, at its most basic, 
the theory attempts to establish convergence of the discrete solution to the true 
solution in an appropriate sense as the discretization parameter tends to zero. 

Underlying the analysis of numerical methods for PDEs is the realization that 
convergence depends on consistency and stability of the discretization. Consistency, 
whose precise definition depends on the particular PDE problem and the type of 
numerical method, aims to capture the idea that the operators and data defining 
the discrete problem are appropriately close to those of the true problem for small 
values of the discretization parameter. The essence of stability is that the discrete 
problem is well-posed, uniformly with respect to the discretization parameter. Like 
well-posedness of the PDE problem itself, stability can be very elusive. One might 
think that well-posedness of the PDE problem, which means invertibility of the 
operator, together with consistency of the discretization, would imply invertibility 
of the discrete operator, since invertible operators between a pair of Banach spaces 
form an open set in the norm topology. But this reasoning is bogus. Consistency is 
not and cannot be defined to mean norm convergence of the discrete operators to the 
PDE operator, since the PDE operator, being an invertible operator between infinite 
dimensional spaces, is not compact and so is not the norm limit of finite dimensional 
operators. In fact, in the first part of the preceding century, a fundamental, and 
initially unexpected, realization was made: that a consistent discretization of a 
well-posed problem need not be stable HHJ 123] ■ Only for very special classes 
of problems and algorithms does well-posedness at the continuous level transfer to 
stability at the discrete level. In other situations, the development and analysis 
of stable, consistent algorithms can be a challenging problem, to which a wide 
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array of mathematical techniques have been applied, just as for establishing the 
well-posedness of PDEs. 

In this paper we will consider PDEs that are related to differential complexes, 
for which de Rham cohomology and Hodge theory are key tools for exploring the 
well-posedness of the continuous problem. These are linear elliptic PDEs, but they 
are a fundamental component of problems arising in many mathematical models, 
including parabolic, hyperbolic, and nonlinear problems. The finite element exterior 
calculus, which we develop here, is a theory which was developed to capture the 
key structures of de Rham cohomology and Hodge theory at the discrete level and 
to relate the discrete and continuous structures, in order to obtain stable finite 
element discretizations. 

1.1. The finite element method. The finite element method, whose develop- 
ment as an approach to the computer solution of PDEs began over 50 years ago 
and is still flourishing today, has proven to be one of the most important tech- 
nologies in numerical PDEs. Finite elements not only provide a methodology to 
develop numerical algorithms for many different problems, but also a mathematical 
framework in which to explore their behavior. They are based on a weak or vari- 
ational form of the PDE problem, and fall into the class of Galerkin methods, in 
which the trial space for the weak formulation is replaced by a finite dimensional 
subspace to obtain the discrete problem. For a finite element method, this subspace 
is a space of piecewise polynomials defined by specifying three things: a simplicial 
decomposition of the domain, and, for each simplex, a space of polynomials called 
the shape functions, and a set of degrees of freedom for the shape functions, i.e., a 
basis for their dual space, with each degree of freedom associated to a face of some 
dimension of the simplex. This allows for efficient implementation of the global 
piecewise polynomial subspace, with the degrees of freedom determining the degree 
of interelement continuity. 

For readers unfamilar with the finite element method, we introduce some basic 
ideas by considering the approximation of the simplest two-point boundary value 
problem 

(1) -u"{x) = fix), -l<a;<l, = =0. 

Weak solutions to this problem are sought in the Sobolev space 1, 1) consisting 
of functions in 1, 1) whose first derivative also belong to L^( — 1, 1). Indeed, 
the solution u can be characterized as the minimizer of the energy functional 

(2) J{u):=^j \u'{x)fdx-j f{x)u{x)dx 

over the space -ff^(— 1, 1) (which consists of 1, 1) functions vanishing at ±1), 

or, equivalently, as the solution of the weak problem: Find u G H^{—1, 1) such that 

(3) J u'{x)v'{x)dx = J f{x)v{x)dx, v € H^{-1,1). 

It is easily seen by integrating by parts that a smooth solution of the boundary value 

problem satisfies the weak formulation and that a solution of the weak formulation 
which posesses appropriate smoothness will be a solution of the boundary value 
problem. 
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Letting Vh denote a finite dimensional subspace of H^{—1,1), called the trial 
space, we may define an approximate solution m/j G as the minimizer of the 
functional J over the trial space (the classical Ritz method), or, equivalently by 
Galerkin's method, in which e Vh is determined by requiring that the variation 
given in the weak problem hold only for functions in Vh, i.e., by the equations 

J u'h{x) v' {x) dx — j f{x)v{x)dx, v€Vh- 

By choosing a basis for the trial space Vh , the Galerkin method reduces to a linear 
system of algebraic equations for the coefficients of the expansion of Uh in terms 
of the basis functions. More specifically, if we write Uh = 'Yl!j=i^j4>ji where the 
functions 4>j form a basis for the trial space, then the Galerkin equations hold if 
and only if Ac — b where the coefficient matrix of the linear system is given by 
Aij = ^^-^4)'j (j)[dx and hi ~ ^^^f<pidx. Since this is a square system of linear 
equations, it is nonsingular if and only if the only solution for f — Q is Uh = 0. This 
follows immediately by choosing v = Uh- 

The simplest finite element method is obtained applying Galerkin's method with 
the trial space Vh consisting of all element of H^{—\, 1) which are linear on each 
subinterval of some chosen mesh — l = a:;o < xi < ••• < xn — lof the domain 
(—1,1). Figure [O] compares the exact solution u — cos{'kx/2) and the finite 
element solution Uh in the case of a uniform mesh with = 14 subintervals. 
The derivatives are compared as well. For this simple problem, Uh is simply the 
orthogonal projection of u into Vh with respect to the inner product defined by the 
left hand side of (|3| , and the finite element method gives good approximation even 
with a fairly coarse mesh. Higher accuracy can easily be obtained by using a finer 
mesh or piecewise polynomials of higher degree. 



Figure 1.1. Approximation of —u" = / by the simplest finite 
element method. The left plot shows u and the right plot shows 
—u', with the exact solution in blue and the finite element solution 
in green. 

The weak formulation ([3| associated to minimization of the functional Q is not 
the only variational formulation that can be used for discretization, and in more 
complicated situations other formulations may bring important advantages. In this 
simple situation, we may, for example, start by writing the differential equation 
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—u" = f as the first order system 

The pair (a, u) can then be characterized variationally as the unique critical point 
of the functional 

I{a,u) — J {-a'^ ~ ua') dx + J fudx 

over H^{—1,1) x 1,1). Equivalently, the pair is the solution of the weak 
formulation: Find a E H^{—1, l),u e 1, 1) satisfying 

J ardx- J ur'dx^O, t e H^i-1,1), 

J a'vdx = J fvdx, uGi^(— 1,1). 

This is called a mixed formulation of the boundary value problem. Note that for the 
mixed formulation, the Dirichlet boundary condition is implied by the formulation, 
as can be seen by integration by parts. Note also that in this case the solution is a 
saddle point of the functional /, not an extremum: /(tr, v) < I{a, u) < I(t, u) for 

T e iJi(-i,i), V e l2(-i,i). 

Although the mixed formulation is perfectly well-posed, it may easily lead to a 
discretization which is not. If we apply Galerkin's method with a simple choice of 
trial subspaces C 1, 1) and C L^{~1, 1), we obtain a finite-dimensional 

linear system, which, however, may be singular, or may become increasingly un- 
stable as the mesh is refined. This concept will be formalized and explored in the 
next section, but the result of such instability is clearly visible in simple compu- 
tations. For example, the choice of continuous piecewise linear functions for both 
E/i and Vh leads to a singular linear system. The choice of continuous piecewise 
linear functions for S/i and piecewise constants for Vh leads to stable discretization 
and good accuracy. However choosing piecewise quadratics for S/j and piecewise 
constants for Vh gives a nonsingular system but unstable approximation (see |25j 
for further discussion of this example). The dramatic difference between the stable 
and unstable methods can be seen in Figure [T^ 

In one dimension, finding stable pairs of finite dimensional subspaces for the 
mixed formulation of the two-point boundary value problem is easy. For any integer 
r > 1, the combination of continuous piecewise polynomials of degree at most r for 
a and arbitrary piecewise polynomials of degree at most r—1 for u is stable as can be 
verified via elementary means (and which can be viewed as a very simple application 
of the theory presented in this paper). In higher dimensions, the problem of finding 
stable combinations of elements is considerably more complicated. This is discussed 



in Section 2.3.1 below. In particular, we shall see that the choice of continuous 
piecewise linear functions for a and piecewise constant functions for u is not stable 
in more than one dimension. However stable element choices are known for this 
problem and again may be viewed as a simple application of the finite element 
exterior calculus developed in this paper. 

1.2. The contents of this paper. The brief introduction to the finite element 
method just given will be continued in Section [2] In particular, there we formalize 
the notions of consistency and stability and establish their relation to convergence. 
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Figure 1.2. Approximation of the mixed formulation for — u" = / 
in one dimension with two choices of elements, piecewise constants 
for u and piecewise linears for a (a stable method, shown in green), 
or piecewise constants for u and piecewise quadratics for cr (unsta- 
ble, shown in red). The left plot shows u and the right plot shows 
tj, with the exact solution in blue. (In the right plot, the blue curve 
essentially coincides with the green curve and hence is not visible.) 

We shall also give several more computational examples. While seemingly simple, 
some of these examples may be surprising even to specialists, and they illustrate 
the difficulty in obtaining good methods and the need for a theoretical framework 
in which to understand such behaviors. 

Like the theory of weak solutions of PDEs, the theory of finite element methods 
is based on functional analysis and takes its most natural form in a Hilbert space 
setting. In Section [3] of this paper, we develop an abstract Hilbert space framework 
which captures key elements of Hodge theory, and can be used to explore the stabil- 
ity of finite element methods. The most basic object in this framework is a cochain 
complex of Hilbert spaces, referred to as a Hilbert complex. Function spaces of such 
complexes will occur in the weak formulations of the PDF problems we consider, 
and the differentials will be differential operators entering into the PDF problem. 
The most canonical example of a Hilbert complex is the de Rham complex of 
a Riemannian manifold, but it is a far more general object with other important 
realizations. For example, it allows the definition of spaces of harmonic forms and 
the proof that they are isomorphic to the cohomology groups. A Hilbert complex 
includes enough structure to define an abstract Hodge Laplacian, defined from a 
variational problem with a saddle point structure. However, for these problems to 
be well-posed, we need the additional property of a closed Hilbert complex, i.e., 
that the range of the differentials are closed. 

In this framework, the finite element spaces used to compute approximate so- 
lutions are represented by finite dimensional subspaces of the spaces in the closed 
Hilbert complex. We identify two key properties of these subspaces: first, they 
should combine to form a subcomplex of the Hilbert complex, and, second, there 
should exist a bounded cochain projection from the Hilbert complex to this sub- 
complex. Under these hypotheses and a minor consistency condition, it is easy to 
show that the subcomplex inherits the cohomology of the true complex, i.e., that 
the cochain projections induce an isomorphism from the space of harmonic forms to 
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the space of discrete harmonic forms, and to get an error estimate on the difference 
between a harmonic form and its discrete counterpart. In the apphcations, this will 
be crucial for stable approximation of the PDEs. In fact, a main theme of finite 
element exterior calculus is that the same two assumptions, the subcomplex prop- 
erty and the existence of a bounded cochain projection, are the natural hypotheses 
to establish the stability of the corresponding discrete Hodge Laplacian, defined by 
the Galerkin method. 

In Section [4] we look in more depth at the canonical example of the de Rham 
complex for a bounded domain in Euclidean space, beginning with a brief summary 
of exterior calculus. We interpret the de Rham complex as a Hilbert complex and 
discuss the PDEs most closely associated with it. This brings us to the topic of 
Section [5] the construction of finite element de Rham subcomplexes, which is the 
heart of finite element exterior calculus and the reason for its name. In this sec- 
tion, we construct finite element spaces of differential forms — piecewise polynomial 
spaces defined via a simplicial decomposition and specification of shape functions 
and degrees of freedom — which combine to form a subcomplex of the de Rham 
complex admitting a bounded cochain projection. First we construct the spaces 
of polynomial differential forms used for shape functions, relying heavily on the 
Koszul complex and its properties, and then we construct the degrees of freedom. 
We next show that the resulting finite element spaces can be efficiently imple- 
mented, have good approximation properties, and can be combined into de Rham 
subcomplexes. Finally, we construct bounded cochain projections, and, having ver- 
ified the hypotheses of the abstract theory, draw conclusions for the finite element 
approximation of the Hodge Laplacian. 

In the final two sections of the paper, we make other applications of the abstract 
framework. In the last section, we study a differential complex we call the elasticity 
complex, which is quite different from the de Rham complex. In particular, one of its 
differentials is a partial differential operator of second order. Via the finite element 
exterior calculus of the elasticity complex, we have obtained the first stable mixed 
finite elements using polynomial shape functions for the equations of elasticity, with 
important applications in solid mechanics. 

1.3. Antecedents and related approaches. We now discuss some of the an- 
tecedents of finite element exterior calculus and some related approaches. While 
the first comprehensive view of finite element exterior calculus, and the first use of 
that phrase, was in the 2006 paper [8], this was certainly not the first intersection 
of finite element theory and Hodge theory. In 1957, Whitney |HS] published his 
complex of Whitney forms, which is, in our terminology, a finite element de Rham 
subcomplex. Whitney's goals were geometric. For example, he used these forms in 
a proof of de Rham's theorem identifying the cohomology of a manifold defined via 
differential forms (de Rham cohomology) with that defined via a triangulation and 
cochains (simplicial cohomology). With the benefit of hindsight, we may view this, 
at least in principle, as an early application of finite elements to reduce the compu- 
tation of a quantity of interest defined via infinite dimensional function spaces and 
operators, to a finite dimensional computation using piecewise polynomials on a 
triangulation. The computed quantities are the Betti numbers of the manifold, i.e., 
the dimensions of the de Rham cohomology spaces. For these integer quantities, 
issues of approximation and convergence do not play much of a role. The situation 
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is different in the 1976 work of Dodziuk |39j and Dodziuk and Patodi [40^, who con- 
sidered the approximation of the Hodge Laplacian on a Riemannian manifold by 
a combinatorial Hodge Laplacian, a sort of finite difference approximation defined 
on cochains with respect to a triangulation. The combinatorial Hodge Laplacian 
was defined in [39 using the Whitney forms, thus realizing the finite difference 
operator as a sort of finite element approximation. A key result in |39l was a proof 
of some convergence properties of the Whitney forms. In jJD] the authors applied 
them to show that the eigenvalues of the combinatorial Hodge Laplacian converge 
to those of the true Hodge Laplacian. This is indeed a finite element convergence 
result, as the authors remark. In 1978, Miiller 71\ further developed this work and 
used it to prove the Ray-Singer conjecture. This conjecture equates a topological 
invariant defined in terms of the Riemannian structure with one defined in terms of 
a triangulation, and was the original goal of ^SSJHO]- (Cheeger [30 gave a different, 
independent proof of the Ray-Singer conjecture at about the same time.) Other 
spaces of finite element differential forms have appeared in geometry as well, espe- 
cially the differential graded algebra of piecewise polynomial forms on a simplicial 
complex introduced by Sullivan [83l [84] . Baker [15] calls these Sullivan- Whitney 
forms, and, in an early paper bringing finite element analysis techniques to bear 
on geometry, gives a numerical analysis of their accuracy for approximating the 
eigenvalues of the Hodge Laplacian. 

Independently of the work of the geometers, during the 1970s and 1980s nu- 
merical analysts and computational engineers reinvented various special cases of 
the Whitney forms and developed new variants of them to use for the solution of 
partial differential equations on two and three-dimensional domains. In this work, 
naturally, implementational issues, rates of convergence, and sharp estimates played 
a more prominent role than in the geometry literature. The pioneering paper of 
Raviart and Thomas [76], presented at a finite element conference in 1975, proposed 
the first stable finite elements for solving the scalar Laplacian in two dimensions 
using the mixed formulation. The mixed formulation involves two unknown fields, 
the scalar- valued solution, and an additional vector- valued variable representing its 
gradient. Raviart and Thomas proposed a family of pairs of finite element spaces, 
one for each polynomial degree. As was realized later, in the lowest degree case the 
space they constructed for the vector-valued variable is just the space of Whitney 
1-forms, while they used piecewise constants, which are Whitney 2-forms, for the 
scalar variable. For higher degrees, their elements are the higher order Whitney 
forms. In three-dimensions, the introduction of Whitney 1- and 2-forms for finite 
element computations and their higher degree analogues was made by Nedelec [72] 
in 1980, while the polynomial mixed elements which can be viewed as Sullivan- 
Whitney forms were introduced as finite elements by Brezzi, Douglas, and Marini 
PB in 1985 in two dimensions, and then by Nedelec [72j in 1986 in three dimensions. 

In 1988 Bossavit, in a paper in the lEE Transactions on Magnetics [H], made 
the connection between Whitney's forms used by geometers and some of the mixed 
finite element spaces that had been proposed for electromagnetics, inspired in part 
by Kotiuga's Ph.D. thesis in electrical engineering [66]. Maxwell's equations are 
naturally formulated in terms of differential forms, and the computational electro- 
magnetics community developed the connection between mixed finite elements and 
Hodge theory in a number of directions. See, in particular, pTl 1571 ITfl [551 [5^ 170] . 
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The methods we derive here are examples of compatible discretization meth- 
ods, which means that at the discrete level they reproduce, rather than merely 
approximate, certain essential structures of the continuous problem. Other ex- 
amples of compatible discretization methods for elliptic PDEs are mimetic finite 
difference methods [16j [27] including covolume methods [74] and the discrete ex- 
terior calculus [38]. In these methods, the fundamental object used to discretize 
a differential fc-form is typically a simplicial cochain, i.e., a number is assigned to 
each fc-dimensional face of the mesh representing the integral of the fc-form over the 
face. This is more of a finite difference, rather than finite element, point of view, 
recalling the early work of Dodziuk on combinatorial Hodge theory. Since the space 
of /c-dimensional simplicial cochains is isomorphic to the space of Whitney fc-forms, 
there is a close relationship between these methods and the simplest methods of 
the finite element exterior calculus. In some simple cases, the methods even coin- 
cide. In contrast to the finite element approach, these cochain-based approaches do 
not naturally generalize to higher order methods. Discretizations of exterior calcu- 
lus and Hodge theory have also been used for purposes other than solving partial 
differential equations. For example, discrete forms which are identical or closely 
related to cochains or the corresponding Whitney forms play an important role in 
geometric modeling, parameterization, and computer graphics. See for example 

[SnUMllMllHT]. 

1.4. Highlights of the finite element exterior calculus. We close this in- 
troduction by highlighting some of the features that are unique or much more 
prominent in the finite element exterior calculus than in earlier works. 

• We work in an abstract Hilbert space setting that captures the fundamental 
structures of the Hodge theory of Riemannian manifolds, but applies more 
generally. In fact, the paper proceeds in two parts, first the abstract theory 
for Hilbert complexes, and then the application to the de Rham complex 
and Hodge theory and other applications. 

• Mixed formulations based on saddle point variational principles play a 
prominent role. In particular, the algorithms we use to approximate the 
Hodge Laplacian are based on a mixed formulation, as is the analysis of the 
algorithms. This is in contrast to the approach in the geometry literature, 
where the underlying variational principle is a minimization principle. In 
the case of the simplest elements, the Whitney elements, the two meth- 
ods are equivalent. That is, the discrete solution obtained by the mixed 
finite element method using Whitney forms, is the same as obtained by 
Dodziuk's combinatorial Laplacian. However the different viewpoint leads 
naturally to different approaches to the analysis. The use of Whitney forms 
for the mixed formulation is obviously a consistent discretization, and the 
key to the analysis is to establish stability (see the next section for the ter- 
minology). However, for the minimization principle, it is unclear whether 
Whitney forms provide a consistent approximation, because they do not 
belong to the domain of the exterior coderivative, and, as remarked in |40j . 
this greatly complicates the analysis. The results we obtain are both more 
easily proven and sharper. 

• Our analysis is based on two main properties of the subspaces used to dis- 
cretize the Hilbert complex. First, they can be formed into subcomplexes, 
which is a key assumption in much of the work we have discussed. Second, 
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there exist a bounded cochain projection from the Hilbert complex to the 
subcomplex. This is a new feature. In previous work, a cochain projection 
often played a major role, but it was not bounded, and the existence of 
bounded cochain projections was not realized. In fact, they exist quite gen- 
erally (see Theorem |3.7| , and we review the construction for the de Rham 
complex in Section |5.5| 

• Since we are interested in actual numerical computations, it is important 
that our spaces be efficiently implementable. This is not true for all piece- 
wise polynomial spaces. As explained in the next section, finite element 
spaces are a class of piecewise polynomial spaces that can be implemented 
efficiently by local computations thanks to the existence of degrees of free- 
dom, and the construction of degrees of freedom and local bases is an 
important part of the finite element exterior calculus. 

• For the same reason, high order piecewise polynomials are of great im- 
portance, and all the constructions and analysis of finite element exterior 
calculus carries through for arbitrary polynomial degree. 

• A prominent aspect of the finite element exterior calculus is the role of 
two families of spaces of polynomial differentials form, Vr^^ and V^A'^, 
where the index r > 1 denotes the polynomial degree and fc > the form 
degree. These are the shape functions for corresponding finite element 
spaces of differential fc-forms which include, as special cases, the Lagrange 
finite element family, and most of the stable finite element spaces that 
have been used to define mixed formulations of the Poisson or Maxwell's 
equations. The space Vi hJ^ is the classical space of Whitney fc-forms. The 
finite element spaces based on VrA^ are the spaces of Sullivan- Whitney 
forms. We show that for each polynomial degree r, there are 2"~^ ways to 
form these spaces in de Rham subcomplexes for a domain in n dimensions. 
The unified treatment of the spaces VrA^ and V^A'^, particularly their 
connections via the Koszul complex, is new to the finite element exterior 
calculus. 

The finite element exterior calculus unifies many of the finite element methods 
that have been developed for solving PDEs arising in fluid and solid mechanics, 
electromagnetics, and other areas. Consequently, the methods developed here have 
been widely implemented and applied in scientific and commercial programs such 
as GetDP [^, FEniCS [44j, EMSolve [35], deal.II gS], Diffpack [M], Getfem-h-h 
[77] . and NGSolve [7H]. We also note that, as part of a recent programming effort 
connected with the FEniCS project, Logg and Mardal [53] have implemented the 
full set of finite element spaces developed in this paper strictly following the finite 
element exterior framework as laid out here and [8]. 



2. Finite element discretizations 

In this section we continue the introduction to the finite element method begun 
above. We move beyond the case of one dimension, and consider not only the for- 
mulation of the method, but also its analysis. To motivate the theory developed 
later in this paper, we present further examples that illustrate how for some prob- 
lems, even rather simple ones, deriving accurate finite element methods is not a 
straightforward process. 
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2.1. Galerkin methods and finite elements. We consider first a simple prob- 
lem, which can be discretized in a straightforward way, namely the Dirichlet prob- 
lem for Poisson's equation in a polyhedral domain C M": 

(4) - Au = / in r?, u = on dO.. 

This is the generalization to n dimensions of the problem (jlj discussed in the 
introduction, and the solution may again be characterized as the minimizer of an 
energy functional analogous to ^ or as the solution of a weak problem analogous to 
This leads to discretization just as for the one-dimensional case, by choosing a 
trial space Vh C H^{^) and defining the approximate solution Uh G Vu by Galerkin's 
method: 

/ gr&d Uh{x) ■ gra,dv{x) dx = I f{x)v{x)dx, vGVh- 
JQ. Jn 
As in one dimension, the simplest finite element method is obtained by using 
the trial space consisting of all piecewise linear functions with respect to a given 
simplicial triangulation of the domain fi, which are continuous and vanish on dfl. 
A key to the efficacy of this finite element method is the existence of a basis for the 
trial space consisting of functions which are locally supported, i.e., vanish on all 



but a small number of the elements of the triangulation. See Figure 2.1 Because 



of this, the coefficient matrix of the linear system is easily computed and is sparse, 
and so the system can be solved efficiently. 



Figure 2.1. A piecewise linear finite element basis function. 



More generally, a finite element method is a Galerkin method for which the trial 
space Vh is a space of piecewise polynomial functions which can be obtained by 
what is called the finite element assembly process. This means that the space can 
be defined by specifying the triangulation and, for each element T G T^, a space 
of polynomial functions on T called the shape functions, and a set of degrees of 
freedom. By degrees of freedom on T, we mean a set of functionals on the space 
of shape functions, which can be assigned values arbitrarily to determine a unique 
shape function. In other words, the degrees of freedom form a basis for the dual 
space of the space of shape functions. In the case of piecewise linear finite elements, 
the shape functions are of course the linear polynomials on T, a space of dimension 
n + 1, and the degrees of freedom are the n + 1 evaluation functionals p i-^ pix), 
where x varies over the vertices of T. For the finite element assembly process, we 
also require that each degree of freedom be associated to a face of some dimension 
of the simplex T. For example, in the case of piecewise linear finite elements, the 
degree of freedom p t-^ p{x) is associated to the vertex x. Given the triangulation. 
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shape functions, and degrees of freedom, the finite element space Vh is defined as 
the set of functions on (possibly multivalued on the element boundaries) whose 
restriction to any T d Th belongs to the given space of shape functions on T, 
and for which the degrees of freedom are single-valued in the sense that when two 
elements share a common face, the corresponding degrees of freedom take on the 
same value. Returning again to the example of piecewise linear functions, Vh is 
the set of functions which are linear polynomials on each element, and which are 
single-valued at the vertices. It is easy to see that this is precisely the space of 
continuous piecewise linear functions, which is a subspace of H^(fl). As another 
example, we could take the shape functions on T to be the polynomials of degree 
at most 2, and take as degrees of freedom the functions p i— > p{x), x a vertex of T, 
and p I— > J^pds, e an edge of T. The resulting assembled finite element space is the 
space of all continuous piecewise quadratics. The finite element assembly process 
insures the existence of a computable locally supported basis, which is important 
for efficient implementation. 

2.2. Consistency, stability, and convergence. We now turn to the important 
problem of analyzing the error in the finite element method. To understand when 
a Galerkin method will produce a good approximation to the true solution, we 
introduce the standard abstract framework. Let F be a Hilbert space, B : V xV ^ 
E a bounded bilinear form, and F : V M. a. bounded linear form. We assume the 
problem to be solved can be stated in the form: Find u Cz V such that 



This problem is called well-posed if for each F e V* , there exists a unique solution 
u € V and the mapping F i—^ u is bounded, or, equivalently, if the operator L : 
V ^ V* given by {Lu,v) — B{u,v) is an isomorphism. For the Dirichlet problem 
for Poisson's equation. 



A generalized Galerkin method for the abstract problem begins with a finite- 
dimensional normed space Vh (not necessarily a subspace oi V), a bilinear form 
Bh ■ Vh X Vh R, and a linear form Fh : Vh ^ M, and defines Uh G Vh by 



A Galerkin method is the special case of a generalized Galerkin method for which Vh 
is a subspace of V and the forms Bh and Fh are simply the restrictions of the forms 
B and F to the subspace. The more general setting is important since it allows 
the incorporation of additional approximations, such as numerical integration to 
evaluate the integrals, and also allows for situations in which Vh is not a subspace 
of V. Although we do not treat approximations such as numerical integration in 
this paper, for the fundamental discretization method we study, namely the mixed 
method for the abstract Hodge Laplacian introduced in Section [Xl} the trial space 
Vh is not a subspace of V, since it involves discrete harmonic forms which will not, 
in general, belong to the space of harmonic forms. 

The generalized Galerkin method ([6| may be written Lh.Uh = Fh where Lh : 
Vh V^ is given by {LhU,v) = Bh{u,v), u,v € Vh- If the finite-dimensional 
problem is nonsingular, then we define the norm of the discrete solution operator, 
||I/^"'^||£(v^- v^), as the stability constant of the method. 



B{u,v)^F{v), veV. 




(6) 



Bh{uh,v) ^ Fh{v), veVh- 
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Of course, in approximating the original problem determined by V, B, and F, 
by the generalized Galerkin method given by Vh, -B/i, and Fii, we intend that the 
space Vh in some sense approximates V and that the discrete forms and 
in some sense approximate B and F. This is the essence of consistency. Our 
goal is to prove that the discrete solution approximates u in an appropriate 
sense (convergence). In order to make these notions precise, we need to compare a 
function in to a function in Vh. To this end, we suppose that there is a restriction 
operator 7r;i : y — > Vh, so that tt/jM is thought to be close to u. Then the consistency 
error is simply LhTThU — Fh and the error in the generalized Galerkin method which 
we wish to control is nhu ~ Uh. We immediately get a relation between the error 
and the consistency error 

TThU -Uh = L^'^iLhT^hU - Fh), 

and so the norm of the error is bounded by the product of the stability constant 
and the norm of the consistency error: 

WiThU-UhWvH < \\Ll'^\\c(v;:,Vu)\\LhT^hU- FhWv;:- 

Stated in terms of the bilinear form Bh, the norm of the consistency error can be 
written 

II r II BhUhU,v) - Fh{v) 
\\LhTThU- FhWv- ^ sup — ^ . 

As for stability, the finite dimensional problem is nonsingular if and only if 

. Bh{u,v) 
7ft mf sup j-— — r—r — > 0, 

o^ueVho^veVn Mv^lMv^ 

and the stability constant is then given by 7^ ^• 

Often we consider a sequence of spaces Vh and forms Bh and Fh where we think 
of /i > as an index accumulating at 0. The corresponding generalized Galerkin 
method is consistent if the Vh norm of the consistency error tends to zero with h 
and it is stable if the stability constant 7,"^ is uniformly bounded. For a consistent, 
stable generalized Galerkin method, \\nhu — Uh\\vh tends to zero, i.e., the method 
is convergent. 

In the special case of a Galerkin method, we can bound the consistency error 

Bh{-KhU,v) - Fh{v) B{nhU-u,v) 

sup = sup -— < \\B\\\\nhU - u\\v. 

a^veVu \\v\\v o^veVh \m\v 

In this case it is natural to choose the restriction tt^, to be the orthogonal projection 
onto Vh, and so the consistency error is bounded by the norm of the bilinear form 
times the error in the best approximation of the solution. Thus we obtain 

WTThU-UhWv <lh^\\B\\ inf \\u-v\\v- 
veVh 

Combining this with the triangle inequality, we obtain the basic error estimate for 
Galerkin methods 

(7) \\u - UhWv < (1 + lh'\\B\\) inf II" - v\\v. 

veVh 

(In fact, in this Hilbert space setting, the quantity in parentheses can be replaced 
with 7^^||i?||, see [HS]-) Note that a Galerkin method is consistent as long as the 
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sequence of subspaces Vh is approximating in V in the sense that 

(8) hni inf \\u - v\\v = 0, u e V. 

h^o veVh 

A consistent, stable Galerkin method converges, and the approximation given by 
the method is quasioptimal, i.e., up to multiphcation by a constant, it is as good 
as the best approximation in the subspace. 

In practice, it can be quite difficult to show that the finite dimensional problem 
is nonsingular and to bound the stability constant, but there is one important case 
in which it is easy. Namely when the form B is coercive, i.e., there is a positive 
constant a for which 

B{v,v) > a\\v\\y, V eV, 

and so > ct- The bilinear form ([s]) for Poisson's equation is coercive, as follows 
from Poincare's inequality. This explains, and can be used to prove, the good 
convergence behavior of the method depicted in Figure 

2.3. Computational examples. 



2.3.1. Mixed formulation of the Laplacian. For an example of a problem that fits 
in the standard abstract framework with a noncoercive bilinear form, we take the 
mixed formulation of the Dirichlet problem for Poisson's equation, already intro- 
duced in one dimension in Section Just as there, we begin by writing Poisson's 
equation as the first order system 

(9) a = — grad u, div cr — f. 

The pair (a, u) can again be characterized variationally as the unique critical point 
(a saddle point) of the functional 

I{(T,u) — f {-a ■ a ~ udiv a) dx + f fudx 
Jn 2 Jfi 

over H{dW-n) x L^{n) where H{div;n) = {a e L^{n) : diva e L^ifl)}. Equiva- 
lently, it solves the weak problem: Find a G H{div; H.), u e L^{il) satisfying 

/ (J ■ T dx ~ I udYVTdx = Q^ T e iJ(div; ri), 
Jn Jn 

/ div av dx ^ / fvdx, v € L^{n). 
Jn Jn 

This mixed formulation of Poisson's equation fits in the abstract framework if 
we define V = H{div; n) x L^{n), 



B{a,u;T,v) = / a ■ t dx — / udivTdx-\- I div crvdx, F{t,v) — / fv 
Jn Jn Jn Jn 



dx. 



In this case the bilinear form B is not coercive, and so the choice of subspaces and 
the analysis is not so simple as for the standard finite element method for Poisson's 
equation, a point we already illustrated in the one-dimensional case. 

Finite element discretizations based on such saddle point variational principles 
are called mixed finite element methods. Thus a mixed finite element for Poisson's 
equation is obtained by choosing subspaces T,h C iJ(div; O) and 14 C L^{Q) and 
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seeking a critical point of / over E/j x Vh- The resulting Galerkin method has the 
form: Find ah G Tih,Uh G Vh satisfying 



ah ■ T dx ~ / Uhdivrdx — O, t ^Y,h, I div ahvdx— / fvdx, v€Vh- 
n Jn Jn Jn 

This again reduces to a linear system of algebraic equations. 

Since the bilinear form is not coercive, it is not automatic that the linear system 
is nonsingular, i.e., that for / = 0, the only solution is ah — 0, Uh — 0- Choosing 
T — ah and v — Uh and adding the discretized variational equations, it follows 
immediately that when f — Q, ah — However, Uh need not vanish unless the 
condition that ^Q^Uhdiv t dx — for all r G implies that Uh = 0. In particular, 
this requires that dim(divS;i) > dimV/j. Thus, even nonsingularity of the approxi- 
mate problem depends on a relationship between the two finite dimensional spaces. 
Even if the linear system is nonsingular, there remains the issue of stability, i.e., a 
uniform bound on the inverse operator. 

As mentioned earlier, the combination of continous piecewise linear elements for 
a and piecewise constants for u is not stable in two dimensions. The simplest stable 
elements use the piecewise constants for u, and the lowest order Raviart-Thomas 
elements for a. These are finite elements defined with respect to a triangular 
mesh by shape functions of the form (a + c + 62:2) and one degree of freedom 
for each edge e, namely a ^ J^a ■ nds. We show in Fig. 2.2 below two numerical 
computations that demonstrate the difference between an unstable and stable choice 
of elements for this problem. The stable method accurately approximates the true 
solution u = x{l — x)y{l — y) on (0, 1) x (0, 1) with a piecewise constant, while the 
unstable method is wildly oscillatory. 




Figure 2.2. Approximation of the mixed formulation for Pois- 
son's equation using piecewise constants for u and for a either 
continuous piecewise linears (left), or Raviart-Thomas elements 
(right). The plotted quantity is u in each case. 



This problem is a special case of the Hodge Laplacian with fc = n as discussed 
briefly in Section 4.2 see especially Section 4.2.4 The error analysis for a variety of 
finite element methods for this problem, including the Raviart-Thomas elements, is 
thus a special case of the general theory of this paper, yielding the error estimates 
in Section [5^ 
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2.3.2. The vector Laplacian on a nonconvex polygon. Given the subtlety of finding 
stable pairs of finite element spaces for the mixed variational formulation of Pois- 
son's equation, we might choose to avoid this formulation, in favor of the standard 
formulation, which leads to a coercive bilinear form. However, while the standard 
formulation is easy to discretize for Poisson's equation, additional issues arise al- 
ready if we try to discretize the vector Poisson equation. For a domain il in 
with unit outward normal n, this is the problem 

(10) — grad div u + curl curl u = /, in fi, u-n — Q^ curlwxn^O, on 917. 

The solution of this problem can again be characterized as the minimizer of an 
appropriate energy functional, 

(11) J{u)^- [ {\dWu\^ + \ curl u\^)dx~ [ f-udx, 

but this time over the space H{curl;il) n iJ(div;il), where _ff (curl; fi) — {u E 
L^{n) I curlii e L^{^)} and (div; f]) = {u e H{div;fl)\u- n = on dfl} with 
H{div; n) defined above. This problem is associated to a coercive bilinear form, 
but a standard finite element method based on a trial subspace of the energy space 
i?(curl; O) niJ(div;J7), e.g., using continuous piecewise linear vector functions, 
is very problematic. In fact, as we shall illustrate shortly, if the domain is a 
nonconvex polyhedron, for almost all f the Galerkin method solution will converge 
to a function that is not the true solution of the problem! The essence of this 
unfortunate situation is that any piecewise polynomial subspace of i?(curl; fl) D 
H{div; Q) is a subspace of H^{Q)nH{div; fi), and this space is a closed subspace of 
if (curl; fi) n H{div; fi). For a nonconvex polyhedron, it is a proper closed subspace 
and generally the true solution will not belong to it, due to a singularity at the 
reentrant corner. Thus the method, while stable, is inconsistent. For more on this 
example, see [35] . 

An accurate approximation of the vector Poisson equation can be obtained from 
a mixed finite element formulation, based on the system: 

(T = — divM, grader + curl curlw = / in il, u-n — 0, curlu x n — on dCl. 

Writing this system in weak form, we obtain the mixed formulation of the problem: 
Find a e H^{fl), u e ff(curl; fi) satisfying 

[ardx- fu.,r.drdx = 0, r e H^n), 
Jn Jn 

/ grad cr • w da; + / curl m • curl v da; = / f-vdx, w G i7(curl; f2). 
Jn Jn Jn 

In contrast to a finite element method based on minimizing the energy (111, a 
finite element approximation based on the mixed formulation uses separate trial 
subspaces of 7? ^ (fi) and (curl; fi) , rather than a single subspace of the intersection 
i?(curl;f])ni?(div;f]). 

We now illustrate the nonconvergence of a Galerkin method based on energy 
minimization and the convergence of one based on the mixed formulation, via com- 
putations in two space dimensions (so now the curl of a vector u is the scalar 
du2/dxi — dui/dx2)- For the trial subspaces we make the simplest choices: for 
the former method we use continuous piecewise linear functions and for the mixed 
method we use continuous piecewise linear functions to approximate a G H^{fl) 
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and a variant of the lowest order Raviart-Thomas elements, for which the shape 
functions are the infinitesimal rigid motions (a — bx2, c + bxi) and the degrees of 
freedom are the tangential moments u t-^ J^u ■ s ds for e an edge. The discrete so- 
lutions obtained by the two methods for the problem when / = (—1,0) are shown 



in Figure 2.3 As we shall show later in this paper, the mixed formulation gives an 
approximation that provably converges to the true solution, while, as can be seen 
from comparing the two plots, the first approximation scheme gives a completely 
different (and therefore inaccurate) result. 





Figure 2.3. Approximation of the vector Laplacian by the stan- 
dard finite element method (left) and a mixed finite element 
method (right). The former method totally misses the singular 
behavior of the solution near the reentrant corner. 



This problem is again a special case of the Hodge Laplacian, now with k = 1. 
See Section |4.2.2[ The error analysis thus falls within the theory of this paper. 



yielding estimates as in Section 5.6 



2.3.3. The vector Laplacian on an annulus. In the example just considered, the 
failure of a standard Galerkin method based on energy minimization to solve the 
vector Poisson equation was related to the reentrant corner of the domain and the 
resulting singular behavior of the solution. A quite different failure mode for this 
method occurs if we take a domain which is smoothly bounded, but not simply 
connected, e.g., an annulus. In that case, as discussed below in Section [331 the 



boundary value problem ( 10 1 is not well-posed except for special values of the forc- 
ing function /. In order to obtain a well-posed problem, the differential equation 
should be solved only modulo the space of harmonic vector fields (or harmonic 1- 
forms) — which is a one-dimensional space for the annulus — and the solution should 
be rendered unique by enforcing orthogonality to the harmonic vector fields. If we 
choose the annulus with radii 1/2 and 1, and forcing function / = (0, x), the result- 
ing solution, which can be computed accurately with a mixed formulation falling 
within the theory of this paper, is displayed on the right in Figure 2.4 However, 



the standard Galerkin method does not capture the non-uniqueness and computes 
the discrete solution shown on the left of the same figure, which is dominated by an 
approximation of the harmonic vector field, and so is nothing like the true solution. 
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Figure 2.4. Approximation of the vector Laplacian on an annu- 
lus. The true solution shown here on the right in an (accurate) 
approximation by a mixed method. It is orthogonal to the har- 
monic fields, and satisfies the difl^erential equation only modulo 
harmonic fields. The standard Galerkin solution using continuous 
piecewise linear vector fields shown on the left, is totally different. 



2.3.4. The Maxwell eigenvalue problem. Another situation where a standard finite 
element method gives unacceptable results, but a mixed method succeeds, arises in 
the approximation of elliptic eigenvalue problems related to the vector Laplacian or 
Maxwell's equation. This will be analyzed in detail later in this paper, and here we 
only present a simple but striking computational example. Consider the eigenvalue 
problem for the vector Laplacian discussed above, which we write in mixed form 
as: Find nonzero (a, u) G H^{il) x _ff(curl; il) and A G M satisfying 



(12) 



fa-rdx- f,r.ar.udx^O, r e H\n). 
Jn Jn 



/ grad a ■ V dx + / curl u ■ curl v dx = X u ■ v dx, v G iJ(curl; f2). 
Jn Jn Jn 

As explained in Section |3.6.1[ this problem can be split into two subproblems. In 

particular, if 7^ m G i?(curl; il) and A G M solves the eigenvalue problem 

(13) / curl w • curl u = A / u-vdx, u G i?(curl; fi), 

Jn Jn 

and A is not equal to zero, then (cr, u), A is an eigenpair for (|12[ ) with a = 0. 

We now consider the solution of the eigenvalue problem ( |13[ ), with two different 

choices of trial subspaces in i?(curl; Q). Again, to make our point it is enough to 



consider a two dimensional case, and we consider the solution of ( 13 1 with f2 a 
square of side length tt. For this domain, the positive eigenvalues can be computed 
by separation of variables. They are of the form + v? with m and n integers: 



1, 1, 2, 4, 4, 5, 5, 8, . . .. If we approximate ( 13 1 using the space of continuous piecewise 



linear vector fields as the trial subspace of i?(curl; fi), the approximation fails badly. 
This is shown for an unstructured mesh in Figure [275| and for a structured crisscross 
mesh in Figure |2.6| where the nonzero discrete eigenvalues are plotted. Note the 
very different mode of failure for the two mesh types. For more discussion of the 
spurious eigenvalues arising using continuous piecewise linear vector fields on a 
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crisscross mesh see [20|. By contrast, if we use the lowest order Raviart-Thomas 



approximation of _ff(curl;il), as shown on the right of Figure 2.5 we obtain a 
provably good approximation for any mesh. This is a very simple case of the 



general eigenvalue approximation theory presented in Section 3.6 below 




Figure 2.5. Approximation of the nonzero eigenvalues of ( 13 1 on 



an unstructured mesh of the square (left) using continuous piece- 
wise linear finite elements (middle) and Raviart-Thomas elements 
(right). For the former, the discrete spectrum looks nothing like 
the true spectrum, while for the later it is very accurate. 








Figure 2.6. Approximation of the nonzero eigenvalues of (13 1 
using continuous piecewise linear elements on the structured mesh 
shown. The first seven discrete nonzero eigenvalues converge to 
true eigenvalues, but the eighth converges to a spurious value. 



3. HiLBERT COMPLEXES AND THEIR APPROXIMATION 

In this section, we construct a Hilbert space framework for finite element exterior 
calculus. The most basic object in this framework is a Hilbert complex, which ex- 
tracts essential features of the de Rham complex. Just as the Hodge Laplacian 
is naturally associated with the de Rham complex, there is a system of variational 
problems, which we call the abstract Hodge Laplacian, associated to any Hilbert 
complex. Using a mixed formulation we prove that these abstract Hodge Laplacian 
problems are well-posed. We next consider approximation of Hilbert complexes 
using finite dimensional subspaces. Our approach emphasizes two key properties, 
the subcomplex property and the existence of bounded cochain projections. These 
same properties prove to be precisely what is needed both to show that the ap- 
proximate Hilbert complex accurately reproduces geometrical quantities associated 
to the complex, like cohomology spaces, and also to obtain error estimates for the 
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approximation of the abstract Hodge Laplacian source and eigenvalue problems, 
which is our main goal in this section. In the following section of the paper we will 
derive finite element subspaces in the concrete case of the de Rham complex and 
verify the hypotheses needed to apply the results of this section. 

Although the de Rham complex is the canonical example of a Hilbert complex, 
there are many others. In this paper, in Section [6] we consider some variations of 
de Rham complex that allow us to treat more general PDEs and boundary value 
problems. In the final section we briefly discuss the equations of elasticity, for which 
a very different Hilbert complex, in which one of the differentials is a second order 
PDE, is needed. Another useful feature of Hilbert complexes, is that a subcomplex 
of Hilbert complex is again such, and so the properties we establish for them apply 
not only at the continuous, but also at the discrete level. 

3.1. Basic definitions. We begin by recalling some basic definitions of homolog- 
ical algebra and functional analysis and establishing some notation. 

3.1.1. Cochain complexes. Consider a cochain complex (V, d) of vector spaces, i.e., 
a sequence of vector spaces V''' and linear maps d*"' , called the differentials: 

. . . y*^^! d''-^ ^ yk yk+l ^ . . . 

with d'^ o d'^^^ — 0. Equivalently, we may think of such a complex as the graded 
vector space V — ® V''" , equipped with a graded linear operator d : V ^ V oi 
degree +1 satisfying d o d = 0. A chain complex differs from a cochain complex 
only in that the indices decrease. All the complexes we consider are nonnegative 
and finite, meaning that V'^ — whenever k is negative or sufficiently large. 

Given a cochain complex {V,d), the elements of the null space 3*^ — 5''{V,d) of 
d'^ are called the fc-cocycles and the elements of the range 03*^ — ^^{V,d) of d^~^ 
the fc-coboundaries. The fcth cohomology group is defined to be the quotient space 

Given two cochain complexes (V, d) and {V , d'), a set of linear maps /''' : V'^ 
V'^ satisfying d'^ = f^'^^d^ (i.e., a graded linear map f : V V oi degree 
satisfying d' f = fd) is called a cochain map. When / is a cochain map, 
maps fc-cocycles to A:-cocycles and fe-coboundaries to fc-coboundaries, and hence 
induces a map / on the cohomology spaces. This map is functorial, i.e., it respects 
composition. 

Let (V, d) be a cochain complex and {Vh, d) a subcomplex. In other words, V/^ is 
a subspace of V'' and d'^V^ C V,^^^. Then the inclusion ih : Vh ^ V is a, cochain 
map and so induces a map of cohomology. If there exists a cochain projection of 
V onto Vfi, i.e., a cochain map iTh such that : V'^ ^ leaves the subspace Vj^ 
invariant, then TTh o i^ — idy^, so ith ° ih = i'^S^/Bj' (where := 'i^{Vh,d) and 
similarly for SfJ. We conclude that in this case i^ is injective and tt^ is surjective. 
In particular, the dimension of the cohomology spaces of the subcomplex is at most 
that of the supercomplex. 

3.1.2. Closed operators on Hilbert space. This material can be found in many places, 
e.g., [Ml Chapter III, § 5 and Chapter IV, § 5.2] or [Ml Chapter II, § 6 and Chapter 
VII]. 

By an operator T from a Hilbert space X to a Hilbert space Y , we mean a linear 
operator from a subspace V oi X, called the domain of T, into Y . The operator 
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T is not necessarily bounded and the domain V is not necessarily closed. We say 
that the operator T is closed if its graph { (a;, Tx) | x € F } is closed in X xY . We 
endow the domain V with the graph norm inner product, 

{v,w)v — {v,w)x + {Tv,Tw)y- 

It is easy to check that this makes V a Hilbert space (i.e., complete), if and only if 
T is closed, and moreover, that T is a bounded operator from V to Y. Of course, 
the null space of T is the set of elements of its domain which it maps to 0, and the 
range of T is T{V). The null space of a closed operator from X to y is a closed 
subspace of X, but its range need not be closed in Y (even if the operator is defined 
on all of X and is bounded). 

The operator T is said to be densely defined if its domain V is dense in X. In 
this case the adjoint operator T* from F to X is defined to be the operator whose 
domain consists of all y Cz Y for which there exists x ^ X with 

{x,v)x ^ {y,Tv)Y, veV, 

in which case T*y = x (well-defined since V is dense). If T is closed and densely 
defined, then T* is as well and T** = T. Moreover, the null space of T* is the 
orthogonal complement of the range of T in Y. Finally, by the closed range theorem, 
the range of T is closed in Y if and only if the range of T* is closed in X. 

If the range of a closed linear operator is of finite codimension, i.e., dim Y/T(V) < 
oo, then the range is closed 60, Lemma 19.1.1]. 

3.1.3. Hilbert complexes. A Hilbert complex is a sequence of Hilbert spaces W'' and 
closed, densely-defined linear operators d'^ from W'^ to W'^+i such that the range 
of d'^ is contained in the domain of d'^^^ and d'''^^ o d^ = 0. A Hilbert complex 
is bounded if, for each A:, d*^ is a bounded linear operator from W'^ to W''^^. In 
other words, a bounded Hilbert complex is a cochain complex in the category of 
Hilbert spaces. A Hilbert complex is closed if for each k, the range of d'' is closed 
in A Fredholm complex is a Hilbert complex for which the range of d'^ 

is finite codimensional in the null space of d'^^^ (and so is closed). Hilbert and 
Fredholm complexes have been discussed by various authors working in geometry 
and topology. Briining and Lesch |28] have advocated for them as an abstraction of 
elliptic complexes on manifolds and applied them to spectral geometry on singular 
spaces. Glotko [55] used them to define a generalization of Sobolev spaces on 
Riemannian manifolds and study their compactness properties and Gromov and 
Shubin [SS] to define topological invariants of manifolds. 

Associated to any Hilbert complex {W,d) is a bounded Hilbert complex {V,d), 
called the domain complex, for which the space V'^ is the domain of d*^, endowed 
with the inner product associated to the graph norm: 

{u,v)vk — {u,v)iYk + {d^U,d*'v)y/k + l. 

Then d^ is a bounded linear operator from to V^^^ , and so (V,d) is indeed a 
bounded Hilbert complex. The domain complex is closed or Fredholm if and only 
if the original complex {W, d) is. 

Of course, for a Hilbert complex (T4^, d), we have the null spaces and ranges 
and 03*^. Utilizing the inner product, we define the space of harmonic forms 
_ 3^ p |.]-^g orthogonal complement of 05 in 3'^- It is isomorphic to 

the reduced cohomology space 3''/*B''' or, for a closed complex, to the cohomology 
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space 'i^ /^'^ . For a closed Hilbert complex, we immediately obtain the Hodge 
decomposition 

(14) W'= = «B'=®io'=®3''^"'- 

For the domain complex (V, d), the null space, range, and harmonic forms are the 
same spaces as for the original complex, and the Hodge decomposition is 

The third summand, 'i^^^ = S'^-L"' n . 

Continuing to use the Hilbert space structure, we define the dual complex {W, d*), 
which is a Hilbert chain complex rather than cochain complex. The dual complex 
uses the same spaces W^, with the differential dl being the adjoint of d''~^, so d*f. is 
a closed, densely-defined operator from W'' to W''~^, whose domain we denote by 
V^. The dual complex is closed or bounded if and only if the original complex is. 
We denote by 31 = S'^-^"' the null space of dl, and by the range of dl_^_^. Thus 
_ ^k Q 2* jg i-jjg space of harmonic forms both for the original complex and the 



dual complex. Since 3*^^"' = 25^, the Hodge decomposition (14 1 can be written 

(15) W'' ^ ® 9)'' ® ^l- 

We henceforth simply write 3'^"'" for 3*^''"^ • 

Let {W,d) be a closed Hilbert complex with domain complex {V,d). Then d'' is 
a bounded bijection from 3'^"'" to the Hilbert space <8'^+^ and hence, by Banach's 
bounded inverse theorem, there exists a constant cp such that 

(16) llt'lk <cp||d'=i;||w, «e3^^, 

which we refer to as a Poincare inequality. We remark that the condition that S'^ 



is closed is not only sufficient to obtain ( 16 1 , but also necessary. 
The subspace V'' D V^* of W'' is a Hilbert space with the norm 

Mlnv = Ml + Mv' = nv\\w + Wd'vllw + Wdlvf^, 

and is continuously included in W'^ . We say that the Hilbert complex ( W, d) has 
the compactness property if V'' H is dense in W'^ and the inclusion is a compact 
operator. Restricted to the space io'^ of harmonic forms, the V'^ n norm is equal 
to the W'^ norm (times \/2). Therefore the compactness property implies that the 
inclusion of ij'' into itself is compact, so Sj'' is finite dimensional. In summary, for 
Hilbert complexes, compactness property Fredholm closed. 

3.2. The abstract Hodge Laplacian and the mixed formulation. Given a 
Hilbert complex {W,d), the operator L = dd* + d*d is an unbounded operator 
W'^ called, in the case of the de Rham complex, the Hodge Laplacian. We 

refer to it as the abstract Hodge Laplacian in the general situation. Its domain is 

= { u e y'^ n Vfc* I e V^+^,d*u e v''-^ }. 

If u solves Lu — /, then 

(17) + (d*M,rf*w) = weF'^nVfc*. 

Note that, in this equation, and henceforth, we use (■ , ■) and || • || without sub- 
scripts, meaning the inner product and norm in the appropriate space. 

The harmonic functions measure the extent to which the Hodge Laplace problem 



(17 1 is well-posed. The solutions to the homogeneous problem (/ — 0) are precisely 
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the functions in Sj'^ . Moreover, a necessary condition for a solution to exist for 
non-zero / e W'^' is that f ± Sj'^. 



For computational purposes, a formulation of the Hodge Laplacian based on (17) 
may be problematic, even when there are no harmonic forms, because it may not be 
possible to construct an efficient finite element approximation for the space V'^OV^ . 
We have already seen an example of this in the discussion of the approximation of 
a boundary value problem for the vector Laplacian in Section |2.3.2[ Instead we 
introduce another formulation, which is a generalization of the mixed formulation 
discussed in Section [2] and which, simultaneously, accounts for the nonuniqueness 
associated with harmonic forms. With {W,d) a Hilbert complex, {V,d) the asso- 
ciated domain complex, and / G W'^ given, we define the mixed formulation of 
the abstract Hodge Laplacian as the problem of finding (cr, u,p) G V^~^ x V'' x Sj'' 
satisfying 

(a,T) - (dr,u) =0, TeV''~\ 
(18) {da,v) + {du,dv) + {v,p) ^ {f,v), ^;ey^ 

{u,q)=0, qeSjK 



Remark. The equations ( 18 1 are the Euler-Lagrange equations associated to a 
variational problem. Namely, if we define the quadratic functional / : V'^^^ xV'' x 
S)'' -> M by 

/(r, v,q) = ^ (r, r) - {dr, - ^ {dv, dv) - {v,q) + (/, v) , 



then a point {a,u,p) E V'^ ^ x V'^ x f)^ is a critical point of / if and only if (18) 
holds, and in this case 

I{cr,v,q) < I{a,u,p) < J{t,u,p), {r.v.q) G V''-^ x V'' x i^^ 

Thus the critical point is a saddle point. 

An important result is that if the Hilbert complex is closed, then the mixed 
formulation is well-posed. The requirement that the complex is closed is crucial, 
since we rely on the Poincare inequality. 

Theorem 3.1. Let {W,d) he a closed Hilbert complex with domain complex {V^d). 
The mixed formulation of the abstract Hodge Laplacian is well-posed. That is, for 
any f € , there exists a unique {(J,u,p) £ V^~^ x V'' x f)^ satisfying ( [T8| . 
Moreover 

\\a\\v + \\u\\v + \\p\\ <c\\f\\, 



where c is a constant depending only on the Poincare constant cp in (16) 



We shall prove Theorem |3.1| in Section |3.2.2[ First, we interpret the mixed 
formulation. 

3.2.1. Interpretation of the mixed formulation. The first equation states that u 
belongs to the domain of d* and d*u = a G V''^^. The second equation similarly 
states that du belongs to the domain of d* and d*du = f — p — da. Thus u belongs 
to the domain Dl of L and solves the abstract Hodge Laplacian equation 

Lu^ f ~p. 

The harmonic form p is simply the orthogonal projection P^f of / onto Sj'^, required 
for existence of a solution. Finally the third equation fixes a particular solution, 
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through the condition u _L . Thus Theorem 3.1 estabhshes that for any f £ W 



rk 



there is a unique u e D]^ such that Lu = f — Ps^f and u ± Sj . We define Kf = u, 
so the solution operator K : is a bounded hnear operator mapping into 

Dl. The solution to the mixed formulation is 



The mixed formulation (18 1 is also intimately connected to the Hodge decompo- 



sition (15). Since da £ B , p G f) , and d*du e 05^, the expression / = da+p+d*du 



is precisely the Hodge decomposition of /. In other words 

P<8 = dd*K, Poj- = d*dK, 

where P<s and Pfg. are the T^^'^-orthogonal projections onto *8'^ and 05^, respec- 
tively. We also note that K commutes with d and d* in the sense that 

dKf = Kdf, / e y^ d*Kg = Kd*g, g £ V^* . 

Indeed, if f G V'" and u = Kf, then u e Dl, which implies that du e V''^^ n 
Also d*M e F''"^ so = / - P?,/ - dd*-u e y'^. This shows that du € Dl. 

Clearly 

Ldu = (dd* + d*d)du = dd* du = d(dd* + d*d)u = dLu = d/, 

and both dit and df are orthogonal to harmonic forms. This establishes that du — 
Kdf, i.e., dKf = Kdf. The second equation is established similarly. 

If we restrict the data / in the abstract Hodge Laplacian to an element of 
or of 05*^, we get two other problems which are also of great use in applications. 
The W problem. If / G «B^, then Kf eW^ satisfies 

d*du^f, ^t_L3^ 

while a — d*u = 0, p = P^f = 0. The solution u can be characterized as the 
unique element of 3'^^ such that 

(19) {du,dv) ^ {f,v), v£-i^^. 

and any solution to this problem is a solution of Lu — f, and so is uniquely 
determined. 

The *B problem. If f e ^B*", then u = Kf G S''' satisfies dd*u = f, while p = Pj,/ = 
0. With a = d*u, the pair {a, u) E V^^^ x S'' is the unique solution of 

(20) {<J,T)-{dT,u)=Q, TEV''-\ {d<7,v) = {f,v), VE^\ 

and any solution to this problem is a solution of Lu = f , a = d*u, and so is uniquely 
determined. 

3.2.2. Well-posedness of the mixed formulation. We now turn to the proof of The- 
orem |3.1[ Let B : X X X M. he a, symmetric bounded bilinear form on a Hilbert 
space X which satisfies the inf-sup condition 

• r Bix.y) 
7 :— mt sup - — - — — > 0. 

Then the problem of finding x £ X such that B{x, y) — F{y) for all y E X is 
well-posed: it has a unique solution x for each F e X* , and \\x\\x < 7^"'^ H-F'llx* 
[T^ . The abstract Hodge Laplacian problem ( 18 1 is of this form, where B : [V'^^^ x 



yk ^ ^/cj ^ ^yk ^ x V'^ X S^''] E denotes the bounded bihnear form 

B{a, u,p; T, V, q) — {a, r) — {dr, u) + {da, v) + {du, dv) + {v,p) — {u, q). 
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and F{t,v,p) = (/,«). 

The following theorem establishes the inf-sup condition, and so implies Theo- 
rem |3TT] 

Theorem 3.2. Let (W,d) be a closed Hilhert complex with domain complex 
There exists a constant 7 > 0, depending only on the constant cp in the Poincare 
inequality (16), such that for any {(T,u,p) E V''^^ x V'' xS)'', there exists {T,v,q) S 

yk-l X l/*: X With 

B{a,u,p;T,v,q) > 7(lk||v + + ||p||)(||r||y + + \\q\\). 

Proof. By the Hodge decomposition, u — u<s + uj-j + u±, where U(q — P^u, u^j = 
Ps^u, and u± — P<s*u. Since utg S S*^, u^b — dp, for some p £ 3^^^^- Since 
duj_ — du, we get using ( 16 1 that 

(21) ||/o||y < cpII^sII, ||u_L|k < cp||(iu||, 

where cp > 1 is the constant in Poincares inequality. Let 
1 



(22) 



-peV 



k-l 



u + da + peV, q = p - Us-, e ^j" . 



From (21 1 and the orthogonality of the Hodge decomposition, we have 
(23) ||t||v + \\v\\v + \\q\\ < C{\\a\\v + \\u\\v + \\p\\). 



We also get, from a simple computation using (21 1 and (22 1, that 
B{a,u,p; T,v,q) 



- 2" " 



- 2" " 



\da\\^ + 
-Wdaf 



\du\\^ + 

-\\duf 

-\\duf 



1 



1 



\P\\ 



(o-, P) 



\us,\ 



c 

1 

2c- 



2 ii^sl 
p 



24' 



\da\\ 



\du\\' 



2cl 



\u<s\ 



24 



2 

\v 



The theorem easily follows from this bound and ( 23 1 . 



□ 



We close this section with two remarks. First we note that in fact Theorem 13.2 



establishes more than the well-posedness of the problem ( 18 1 stated in Theorem |3.1| 
It establishes that, for any G € (V''-^)* , F € {V'')* , and R S (^'^)* (these are 
the dual spaces furnished with the dual norms), there exists a unique (cr, u,p) G 
V'^~^ X V'' X Sj'' satisfying 

{a,T)-{dr,u)=G{T), TeV^-\ 
{da,v) + {du,dv) + {v,p) ^ F{v), v e V'' , 

{u,q) = R{q), gei^^ 



and moreover the correspondence (cr, u,p) <-> {F, G, R) is an isomorphism of V 
yk ^ ^k jjjg^j space. 



fc-1 
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Second, we note that the above result bears some relation to a fundamental 
result in the theory of mixed finite element methods, due to Brezzi 24J, which we 
state here. 

Theorem 3.3. Let X and Y be Hilhert spaces and a:XxX^M, 6:XxF^M 
hounded bilinear forms. Let Z ~ {x € X \ b{x,y) — Cz Y }, and suppose that 
there exists positive constants a and 7 such that 

(1) (coercivity in the kernel) a{z,z) > a||z||^, z ^ Z, 

b(x J 

(2) (inf-sup condition) inf sup - — - — — - — > 7. 

^T^yi^yo^xex \\x\\x\\y\\Y 

Then, for all G G X* , F (zY* , there exists a unique u €z X , v Y such that 

a{u, x) + b{x, v) = G{x), x E X, 
^^'^^ b{u,y)^F{y), yeY. 

Moreover \\u\\x + \\v\\y < c(||G|jx' + ll^llr*) with the constant c depending only 
on a, 7, and the norms of the bilinear forms a and b. 

If we make the additional assumption (usually satisfied in applications of this 
theorem), that the bilinear form a is symmetric and satisfies a{x,x) > for all 
^ X G X, then this theorem can be viewed as a special case of Theorem |3.2[ 
In fact, we define W'^ as the completion of X in the inner product given by a, let 
= Y, and define d as the closed linear operator from to with domain 
X given by 

{dx,y)Y = b{x,y), xeX,yeY. 
In this way we obtain a Hilbert complex (with just two spaces W'^ and W^). The 



inf-sup condition of Theorem 3.3 implies that d has closed range, so it is a closed 



Hilbert complex, and so Poincarc's inequality holds. The associated abstract Hodge 



Laplacian is just the system (|24|, and from Theorem 3.2 and the first remark above, 
we have 

.ya{u,u) + \\du\\Y + \\v\\y < c{\\G\\x' + \\F\\y'). 
But, using the coercivity in the kernel, the decomposition X = Z + Z^, and 
Poincare's inequality, we get < C[-y/ a{u, u) + \\du\\Y], which gives the estimate 

from Brezzi's theorem. Finally, we mention that we could dispense with the extra 
assumption about the symmetry and positivity of the bilinear form a, but this 
would require a slight generalization of Theorem |3 . 2| which we do not consider here. 



3.3. Approximation of Hilbert complexes. The remainder of this section will 
be devoted to approximation of quantities associated to a Hilbert complex, such as 
the cohomology spaces, harmonic forms, and solutions to the Hodge Laplacian, by 
quantities associated to a subcomplex. 

Let {W,d) be a Hilbert complex with domain complex {V,d), and suppose we 
choose a finite-dimensional subspace of V'^ for each k. We assume that dVj^ C 
V^'*'"^ so that (Vhjd) is a subcomplex of {V,d). We also take Wj^ to be the same 
subspace but endowed with the norm of W''. In this way we obtain a closed 
(even bounded) Hilbert complex {Wh,d) with domain complex (Vhjd), and all the 
results of Sections |3 . 1 . 3| and |3 .2| apply to this subcomplex. Although the differential 
for the subcomplex is just the restriction of d, and so does not need a new notation, 
its adjoint dl : V^+'^ V^, defined by 

{dlu,v)^{u,dv), ueV^+\veV^, 
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is not the restriction of d* . We use the notations *B/i, 3/i, — ih, ^h, Kh, 
with the obvious meanings. We use the term discrete when we wish to emphasize 
quantities associated to the subcomplex. For example, we refer to f)^ as the space 
of discrete harmonic fc- forms, and the discrete Hodge decomposition is 

The VF^-'-projections P^^ : *B^, P^j : satisfy Pgj^ = ddlKu 

and Prs* = d*f^dKh, respectively, when restricted to Vf^ . We also have that Kh 
commutes with both d and d^ . Note that 05 1 C 25 and 3/5 C 3'^, but in general 
S)^ is not contained in , nor is contained in 3*^^- 

In order that the subspaces Vj^ can be used effectively to approximate quantities 
associated to the original complex, we require not only that they form a subcomplex, 
but of course need to know something about the approximation of by V^*^, i.e., 
an assumption that infj,gy^fc |jit — w|| y is sufficiently small for some or all u S . A 
third assumption, which plays an essential role in our analysis, is that there exists 
a hounded cochain projection from the complex (V, d) to the subcomplex (Vh, d). 
Explicitly, for each k, tt^ maps to Vj^ , leaves the subspace invariant, satisfies 
d^Ti\ = Tr'^^^d'', and there exists a constant c such that ||7r^u||y < cllwllv/ for all 
V e V''. In other words, we have the following commuting diagram relating the 
complex {V,d) to the subcomplex {Vh,d): 

0^V° — ^ — > — - — > ■ ■ ■ — ^ — > ^ 

Note that a bounded projection gives quasioptimal approximation: 

\\u-niiu\\v= inf ||(/ — 7r/i)(M — t;)||y < c inf — 

We now present two results indicating that, under these assumptions, the space 
Sj^ of discrete harmonic forms provides a faithful approximation of Sj'^. In the first 
result we show that a bounded cochain projection into a subcomplex of a bounded 
closed Hilbert complex which satisfies a rather weak approximability assumption 



(namely (25 1 below), induces, not only a surjection, but an isomorphism on coho- 
mology. 

Theorem 3.4. Let {V,d) he a hounded closed Hilhert complex, {Vind) a Hilhert 
suhcomplex, and tt/j a hounded cochain projection. Suppose that for all k 

(25) ||q-7r^<7||y < ||q||y, ^ q ^ S)\ 

Then the induced map on cohomology is an isomorphism. 

Proof. We already know that the induced map is a surjection, so it is sufficient to 
prove that it is an injection. Thus, given z G 3*^ with tt^z G 05^, we must prove 
that z e 25'°. By the Hodge decomposition, z = q + b with q G 9)^ and 6 £ S'^. We 
have that HhZ £ 25^ by assumption and nhb G 25^ since b G 05'° and tt/j is a cochain 



map. Thus n^q = tt/jZ — nhb e 25^^ C 25'°, and so iThq -L q. In view of ( |25| , this 
implies that q — 0, and so z e as desired. □ 
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Remark. In applications, the space of harmonic forms, i^'^, is a finite-dimensional 
space of smooth functions, and tt^ is a projection operator associated to a triangu- 



lation with mesh size h. The estimate (25) will then be satisfied for h sufficiently 



small. However in the most important application, in which (V, d) is the de Rham 
complex and (Vh, d) is a finite element discretization, TTh induces an isomorphism on 
cohomology not only for h sufficiently small, but in fact for all h. See Section [53] 



The second result relating i^'^ and is quantitative in nature, bounding the dis- 
tance, or gap, between these two spaces. Recall that the gap between two subspaces 
E and of a Hilbert space V is defined Chapter IV, § 2.1] 

(26) gap(i?, F) = maxi sup inf — sup inf — 

\\u\\ = l \\v\\ = l 

Theorem 3.5. Let (V,d) be a bounded closed Hilbert complex, (yh,d) a Hilbert 
subcomplex, and tth a bounded cochain projection. Then 

(27) \\{I-P^,)q\\v<\\{I-7:t)q\\v, 9 G ^^ 

(28) \\{I-P^)q\\v<\\{I~T:l)P^q\\v, q ^ 

(29) gap(i3^i3^) < sup \\{I-Trfi)q\\v. 

Il9ll=l 

Proof. Given q G S^^ , Pf,^q = P^^q, since 3^ = and q ± ^'^ D 05^ Also 
TT^q e 3^, since tt/j is a cochain map. This implies (27l. 



If (J € Sj'H C 3^ C 3*^, the Hodge decomposition gives us q — P^q E OS'"', so 
{q — Pi^q) G ^fi, and so is orthogonal to both the discrete harmonic form q and 
3 harmonic form Pf,q. Therefore 

Ik - PfiQlW < Ik - P^Q - ~ P9,q)\\v = \\il- Trh)P9,q\\v- 



Finally (29 1 is an immediate consequence of (27) and (28). □ 



Next we deduce another important property of a Hilbert subcomplex with a 
bounded cochain projection. Since the complex {V,d) is closed and bounded, the 



Poincare inequality (16) holds (with the W and V norms coinciding). Now we 
obtain the Poincare inequality for the subcomplex with a constant that depends 
only on the Poincare constant for the supercomplex and the norm of the cochain 
projection. In the applications, we will have a sequence of such subcomplexes 
related to a decreasing mesh size parameter and this theorem will imply that the 
discrete Poincare inequality is uniform with respect to the mesh parameter, an 
essential step in proving stability for numerical methods. 

Theorem 3.6. Let {V,d) be a bounded closed Hilbert complex, (V/i,rf) a Hilbert 
subcomplex, and tt^ a bounded cochain projection. Then 

\\v\\v <Cp\\TT^\\\\dv\\v, V£5h^, 

where cp is the constant appearing in the Poincare inequality (16) and ||7r^|| denotes 
the V*' operator norm of tt^. 



Proof. Given v S 'i^^, define z g 3'^"'" C by dz = dv. By (16), ||z|| < cp||(iw||, so 
it is enough to show that \\v\\v < ||7r/i|| ||z|| Now, v — tt/^z e V^&nd d{v'-iihz) — Q, 
so V — TTfiZ G ^h- Therefore 

ll^lly = {v,Trhz)v + {v,v~ ■nhz)v = {v,T^hz)v < WAvW'^hAW , 
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and the result follows. □ 

We have estabhshed several important properties possessed by a subcomplex of a 
bounded closed Hilbert complex with bounded cochain projection. We also remark 
that from ( 28 1 and the triangle inequality, we have 

(30) \\q\\v<c\\Pm\\v, qe^l 

We close this section by presenting a converse result. Namely we show that if 
the discrete Poincare inequality and the bound (30) hold, then a bounded cochain 
projection exists. 

Theorem 3.7. Let (V,d) be a bounded closed Hilbert complex and (Vhjd) a sub- 
complex. Assume that 

\\v\\v <ci\\dv\\v, ve'5l-^, and \\q\\v < C2\\Pf,q\\v , q 

for some constants c\ and Ci. Then there exists a bounded cochain projection tt^ 
from {V,d) to (V/j,d), and the V operator norm |j7r/j|| can be bounded in terms of 
ci and C2- 

Proof. As a first step of the proof we define an operator Qh ■ — > 3/^^ by 
dQhV = -PsfeC?'^'- By the first assumption, this operator is T^-bounded since 

IIQht^lly < ci||P>8^(ii;|| < ci||t;||y, 

and liv &Vh then QhV = P^i^v. By the second assumption, the operator Pf, jf,^ has 
a bounded inverse Rh mapping i^'^ to ijj^. Note that for any v G 3fi, RhP^v = Ps^^v- 
We now define nh : V'^ ~* V,f by 

TTh = P^, +RhPf,{i-Qh) + Qh- 

This operator is bounded in V'^ and invariant on V/^, since the three terms cor- 
respond exactly to the discrete Hodge decomposition in this case. Furthermore, 
TThdv = P<Sf^dv = dQhV = diThV, so TTh is indeed a bounded cochain projection. □ 

3.4. Stability and convergence of the mixed method. Next we consider a 
closed Hilbert complex {W, d) and the Galerkin discretization of its Hodge Lapla- 
cian using finite dimensional subspaces of the domain spaces V''. Our main 



assumptions are those of Section 3.3 first, that dV^ C V^'^^, so that we obtain a 
subcomplex 

_> yo ^ y,i ^ . . . ^ v;; ^ 0, 

and, second, that there exists a bounded cochain projection tt/j from (V, d) to (14,, d). 

Let / € W'^. In view of the mixed formulation ( 18 1, we take as an approximation 
scheme: Find ah G V^^^, Uh G V/^ , ph G -ft^, such that 

{ah, t) - (dr, Uh) ==0, re Vl"-^ , 

(31) {dah,v) + {duh,dv) + {v,ph) ^ {f,v), vGVh, 

{uh,q) = 0, qGS)t 

(Recall that we use (•, •) and || • || without subscripts for the W inner product 
and norm.) The discretization (31) is a generalized Galerkin method as discussed 
in Section 2.2 In the case that there are no harmonic forms (and therefore no 
discrete harmonic forms), it is a Galerkin method, but in general not, since S)^ is 
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not in general a subspace of Sj'' . We may write the solution of (31 1, which always 
exists and is unique in view of the results of Section |3.2.2[ as 

where Ph ■ Vj^ is the H^'^-orthogonal projection. 

As in Section 2.2 we will bound the error in terms of the stability of the dis- 
cretization and the consistency error. We start by establishing a lower bound on 
the inf-sup constant, i.e., an upper bound on the stability constant. 

Theorem 3.8. Let (Vh, d) be a family of subcomplexes of the domain complex (V, d) 
of a closed Hilbert complex, parametrized by h and admitting uniformly V -bounded 
cochain projections. Then there exists a constant 7/j > 0, depending only on cp and 
the norm of the projection operators TTh, such thatforany {a,u,p) G ^/f ~^ x X-^^; 



there exists {T,v,q) G ^ >^ ^ '^^^^ 



k-l 

B{a,u,p;T,v,q) > JhihWv + \\u\\v + \\p\\){\\t\\v + \\v\\v + \\q\\)- 



Proof. This is just Theorem 3.2 applied to the Hilbert complex {Vh,d), combined 



with the fact that the constant in the Poincare inequality for is cp ||7r/i|| by 
Theorem [3J1 □ 



From this stability result, we obtain the following error estimate. 

Theorem 3.9. Let (Vh, d) be a family of subcomplexes of the domain complex (V, d) 
of a closed Hilbert complex, parametrized by h and admitting uniformly V -bounded 
cochain projections, and let {<T,u,p) G V^^^ x x ij*^ be the solution of problem 



(I8I and {(Jh,Uh,Ph) G V,f ^ x Vj'^ x the solution of problem (311. Then 

\W ~ (^h\\v + \\u - UhWv + Hp - Phil 

<C{ inf ||(7 — T||y+ inf ||u — w||y+ inf \\p — q\\v + inf \\P<su — v\\v), 
rGV*"^ vevl; qev^ vev^ 



where fi ^ fi'll^ = sup ||(/ - 7rj;)r|| . 
M=i 

Proof. First observe that {a,u,p) satisfies 

B{a,u,p;Th,Vh,qh) = (/, w^) - {Th,Vh,qh) e Vj^'^ x V,f x io^. 

Let r, w, and q be the F-orthogonal projections of tr, u, and p into V^'^^^, , and 
respectively. Then, for any {Th,Vh,qh) G V^~^ x V^' x ij^, we have 

B{ah -T,Uh- v,ph - q;Th,Vh,qh) 

= B{(7 -T,u-v,p- q; Th,Vh,qh) + {u,qh) 

= B{a -T,u-v,p~ q; Th,Vh,qh) + {Pf,^u,qh) 

< Ci\\a -t\\v + \\u - v\\y + \\p-q\\ + ||Pfl,w||)(||T,||v + WvkWv + khW). 
Theorem |3.8| then gives 

(32) ||a-,j -t\\v + \\uh -v\\v + \\ph - q\\ 

< Ci\\a ~ t||v + ||m - ^;||y + lb - 9ll + \\P^,M\)- 



Using (27 1 and the boundedness of the projection we have 



(33) ||p-g|| < ||(/-^,OpII <C inf \\p-q\\v- 

qev^ 
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Next we show that 

Now u _L Sj'^, so u — u<s + uj_, with Ufg G S'^ and uj_ G 3'^"'~- Since C 3'^, 
Pfjh'^± = 0, and since TThU<B £ 23^, Ps,„T^hU<B = 0. Let g = Ps5,^it/||Pf,^u|| G Sj^. 



By Theorem 3.5 there exists r ^ Sj'^ (and so r ± S'') with 

Ik-r|| < ||(/-7r;,)r|l < sup |1 (/ - 7r;,)r|| . 



< 1 and 



llr||=l 



Therefore 

(34) \\Pf,^u\\ = {u^-7ThU<s,q-r) 



< Wil - T^h)u's\\ sup 7r/i)r|l < inf |lP«8W-w|!y) 



lkll=i 



since tt/j is a bounded projection. The theorem follows from (32h(34| and the 
triangle inequality. □ 

To implement the discrete problem, we need to be able to compute the discrete 
harmonic forms. The following lemma shows one way to do this; namely it shows 
that the discrete harmonic forms can be computed as the elements of the null space 
of a matrix. For finite element approximations of the de Rham sequence, which 
are the most canonical example of this theory and which will be discussed below, 
it is often possible to compute the discrete harmonic forms more directly. See, for 
example, [2]. 



Lemma 3.10. 

such that 



Consider the homogeneous linear system: Find (ah, ^h) 



T G V, 



{(Th,T) = {dT,Uh 

{dah,v) + {duh,dv) =0, t; G V/^ 
Then (cr;j,u^) is a solution if and only if Uh ~ 
Proof. Clearly (0,u/i) is a solution if G 



and Ufi G ^^j. 



On the other hand, if {(Th,Uh) is a 
solution, by taking t = au, v = Uh, and combining the two equations, we find that 
||(T/i|p + ||c?u/i||^ = 0, so that (T/i = and duh = 0. Then the first equation implies 
that {dr, Uh) = for all t G V^^^^ , so indeed Uh E Sj'j^. □ 

3.5. Improved error estimates. Suppose we have a family of subcomplexes 
{Vh,d) of the domain complex (V, d) of a closed Hilbert complex, parametrized by h 
with uniformly bounded cochain projections iTh- Assuming also that the subspaces 
are approximating in in the sense of ([s]), we can conclude from Theorem 3.9 
that ah a, Uh ^ u, and ph p a,s h ^ (in the norms of V''~^ and V''). In 
other words, the Galerkin method for the Hodge Laplacian is convergent. 

The rate of convergence will depend on the approximation properties of the 
subspaces V^, the particular component considered {ah, Uh, or ph), the norm in 
which we measure the error (e.g. , W or V) , as well as properties of the data / and the 
corresponding solution. For example, in Section[5]we will consider approximation of 
the de Rham complex using various subcomplexes for which the spaces consist 
of piecewise polynomial differential forms with respect to a triangulation 7^ of the 
domain with mesh size h. The space W'^ is the space of differential fc-forms 
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in this case. One possibility we consider for the solution of the Hodge Laplacian 
for /c-forms using the mixed formulation is to take subspa ces 7 ^r+iA''~^(7/J and 
Vr^^{Th)- The space Vr^^{Th), which is defined in Section 5.2 consists piecewise 
of all fe-forms of polynomial degree at most r. Assuming that the solution u to the 
Hodge Laplacian is sufficiently smooth, application of Theorem 3.9 will give, in this 
case, 

\\(^ ~ <Jh\\v + \\u - Uh\\v + \\p-ph\\ = 0(/i''). 
Approximation theory tells us that this rate is the best possible for ||u — but 
we might hope for a faster rate for \\u — Uh\\ and for \\a — OhWv and \\a — ah\\- 

In order to obtain improved error estimates, we make two addition assumptions, 
first that the complex {yV, d) satisfies the compactness property introduced at the 
end of Section [3.1 1 and second, that the cochain projection is bounded not only in 
V but in W: 

• The intersection n is a dense subset of with compact inclusion. 

• The cochain projections tt^ are bounded in C{W^ ,W^) uniformly with 
respect to h. 

The second property implies that tt/j extends to a bounded linear operator 

. Since the subspaces V/^ are approximating in W'^ as well as in V'' (by density) , 

it follows that tt/j converges pointwise to the identity in W''. Finally, note that if 

we are given a VF-bounded cochain projection mapping W'' — > V/^, the restrictions 

to V'' define a l^-bounded cochain projection. 

Next, we note that on the Hilbert space V'^ n V^T, the inner product given by 

{u,v)vnv' ■— {d*u,d*v) + {du,dv) + {P^u,P^jv) 

is equivalent to the usual intersection inner product which is the sum of the inner 
products for and V^. This can be seen by Hodge decomposing u as P<su + 
PfjU + P<B'U, and using the Poincare bound llPrstull < c||dPf8*'fi|| — c\\du\\ and the 
analogous bound HPswH < c||d*Pf8M|| = c||(i*u||. Now K maps W'^ boundedly into 
V'' n T4* and satisfies 

{Kf,v)vnv* ={f,v-P^v), I^W\ v^V''C^V:. 

In other words, K G C{W^ , n V^) is the adjoint of the operator (/ — Pg^)T G 
C{V'^ n V^,W'') where I G C{V'^ n ,W^) is the compact inclusion operator. 
Hence X is a compact operator W'^ n and, a fortiori, compact as an 

operator from to itself. As an operator on , K is also self-adjoint, since 

(/, Kg) = (f- P^f, Kg) = {dd*Kf + d*dKf, Kg) = (d*Kf, d*Kg) + {dKf, dKg) 

for all /, g G W''. Furthermore, if we follow K by one of the bounded operators 
d : V'' ^ p^fe+i or d* : W'^^^, the compositions dK and d*K are also 

compact operators from W'' to itself. Since we have assumed the compactness 
property, dimi^*"' < oo, and so P^^k is also a compact operator on W'^. Define 

V ^ Vh^ ma^ilK-^ - T^h)dK\\c{w'"i .,w''-^+'-)^ \\(^ ~ '^h)d* K\\c(wi'+^ ,w''+^-^)]- 



(Note that fi already appeared in Theorem 3.9 ) Recalling that composition on the 
right with a compact operator converts pointwise convergence to norm convergence, 
we see that rj, 6, fj, ^ as h ^ 0. In the applications in Sections |5] and |6] the 
spaces A^ will consist of piecewise polynomials. We will then have rj — 0{h), 
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S = 0(/i"""(^'''+^'), and ^ — 0{h^^^) where r denotes the largest degree of complete 
polynomials in the space A{^. 

In Theorem |3.9[ the error estimates were given in terms of the best approximation 
error afforded by the subspaces in the V norm. The improved error estimates will 
be in terms of the best approximation error in the W norm, for which we introduce 
the notation 

E{w)^E^{w)^ inf \\w~v\\, weW''. 

The following theorem gives improved the improved error estimates. Its proof 
incorporates a variety of techniques developed in the numerical analysis literature 
in the last decades, e.g., [4811411 18]. 

Theorem 3.11. Let (V, d) he the domain complex of a closed Hilbert complex 
(W^d) satisfying the compactness property, and {Vh,d) a family of subcomplexes 
parametrized by h and admitting uniformly W -bounded cochain projections. Let 
{(j,u,p) £ V^^^ X X i^*^ be the solution of problem (18) and {(Th,Uh,Ph) G 
V^^^ X V]l^ X the solution of problem (31 1. Then for some constant C indepen- 
dent of h and {a,u,p), we have 

(35) \\d{a-ah)\\<CE{dc7), 

(36) \\a - ahW < C[E{a) + 7]E{da)], 

(37) \\p - PhW < C[Eip) + fiEida)]. 

(38) \\d{u - ui,)\\ < C{E{du) + Tj[Eida) + E{p)]), 
(39) 

\\u - ui,\\ < C{E{u) + 7j[E{du) + E{a)] + (rj^ + S)[E{da) + E{p)] + fiE{P<su)]). 

We now develop the proof of Theorem |3.11| in a series of lemmas. 

Lemma 3.12. Let Vh E 3^^ and v = P<s*Vh- Then 

\\v-Vh\\ < ||(/-7r^)w|| < ri\\dvh\\. 

Proof. Since HhV—Vh G 3^ C 3'^, Vh^v _L TThV~Vh, so, by the Pythagorean theorem, 
— w|| < II (/ — 7r;i)i;|| . The second inequality holds since v = d*Kdvii. □ 



Lemma 3.13. The estimates ( |35[ ) and (36 1 hold. Moreover 
(40) \\P<s,iu-Uh)\\ <C[r^E{a) + ir^'' + 6)Eida)]. 

Proof. Since da^ = P<Shf — P'Sh^'sf — P'Sh'^^^ have 

\\d{a-aH)\\ - ||(/-P<8jda|| < \\{I-TTH)da\\ <CE{da), 



giving (351. To prove ( 36 1 we write 

a = d*Kd(j = d*K{I - P<sjda + d*KP^J(j =: cr^ + cr^ 



Taking t = vr/i cr — cr/j in (18) and (31 1, we obtain 

(cr - (Jh,TTk.a2 - cTh) = (rf(7r/icr2 -ah),u- Uh) = 0. 

Hence, 

h - ^h\\ < Ik - ^h^^W <\\il- ^hM + hi^a'W < CEia)\ + C\\a'\\. 
Since cr^ € S^, 

\\a'f = {Kda\da') = {Kda\ (I - P^,)da) = ((/ - P^jKda\ (/ - P^Jda) 
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< ~ TTh)Kda^\\\\{I - ni,)da\\ < C7j\\<j^\\Eida). 
Then (36 1 follows from the last two estimates. 



Let e — P<Si^(u — Ufi). To estimate e, we set w — Ke, cp = d*w, — K^e, 
(f>h — d'^Wh- Then d(f> = dnhip = d(f>h — e, so tt/i^ — (j)h ^ ^^nd so is orthogonal 
io(l)-^H- Thus U-(l)h\\ < = \\{I-TTh)d*Ke\\. Then 

||e||^ = {d(f>h,e) = {d(f)h,u- Uh) = (o" - (ihAh) 
= (cr - CT/i, (j)h-4>) + {d{<T - (Jh),w) 

(fT - (Th,(j}h -'/'> + {{I - P<Bh)dcr, {I -~ P<Sh)w) 
<\\(7- ahWUh - cl)\\ + 7r,,)dCT||||(/ - 7r,,)u;|| 

< \\a - ahWWil ~ 7:h)d* Ke\\ + CE{da)\\{I - TTh)Ke\\ 

< [7^\\a~ah\\+CE{da)6]\\e\\. 



Combining with ( 36 ) we get ( 40 1 



□ 



Lemma 3.14. Estimate (37l holds and, moreover, |1P^^m|| < Cfj,E{PtQu). 



Proof. The second estimate is just (34 1. Using the Hodge decomposition of / and 
the fact that Pq^S^ = 0, we have Pf,,^p ~ Ph = Pf,uPfi f - P^uf = Paul's-, where 
f<z = P'sf- Therefore, 

\\p-M\ = \\p-P^M\ + \\P^M<CE{p) + \\p^,M\, 



by (271. Applying (28 1 we get 



Il-Ps5h/»II^ = (/S - 7r;,./<B,Pf,„/2}) = (/s - Trhf'B,Pf}Hf'B " PfjPfjnM 

<CUl-7rh).M\\{I-^h)PsiPs^,M<CE{da)fx\\Ps,J^\ 
Combining these results completes the proof of the lemma. 



□ 



Lemma 3.15. The estimate (38) holds 



Proof. From ([18| and (jSlJ, 

(41) {d{a - ah),Vh) + {d{u - Uh),dvh) + {vh,p- Ph) = Vh<EVh. 
Choose Vh = PfB* (T^hU — Uh) in (41 1 and set v ~ d* Kdvh — P<s*Vh. Using Lem- 



mas 



3.12 and (p7|, we get 



{d{u - Uh),d{nhU - Uh)) = (d{u - Uh),dvh) = -(d{c7 - (7,0 + {p-Ph),Vh) 
= -{d{a - ah) + {p- Pfj,^p),Vh ~v) < [\\d{cr - (T,,)|| + \\p - Pi,^p\\]\\vh - v\\ 
< C[E{da) + E{p)]f]\\dvh\\ < Cr][E{da) + E{p)]\\d{7rhU - Uh)\\. 

Hence, 

\\d{TThU - = {d{TrhU ~ u),d{-KhU - Uh)) + {d{u - Uh),d{TThU - Uh)) 

< {\\d{TThU - u)\\ + Cr][E{da) + E{p)]}\\d{TThU - Uh)\\ 

< C{E{du) + v[E{da) + E{p)]}\\d{7:hU - 

The result follows by the triangle inequality. 

Lemma 3.16. \\P<s.^{u - Uh)\\ < C{E{u) + riE{du) + {r]^ + S)[E{d(7) + E{p)])}. 



□ 
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Proof. Again letting Vh = P'^sii'^hU — ut) and v = Pt^'^h, we observe that 

\\P^l{u-Uh)\\ < \\P<si{u-nhu)\\ + \\vh\\ < \\u ~ TThuW + \\vh\\ <CE{u) + \\vh\\. 
Next 

\\Vh\\^ = {vh - V,Vh) + {VjlThU - U,i) = {vh - V,Vh) + {v.TThU ~ u) + {v , U - Uh) ■ 

Then using Lemmas |3.12| and |3.15| we get 

{vh - v,Vh) + {v,TThU - u) < \\vh - v\\\\vh\\ + \\v\\\\{I - 7rft)u|| 

< [vWdVhW + ~ TTh)u\\]\\vh\\ < [Tj\\d{TThU ~ Uh)\\ + CE{u)]\\vh\\. 

We next estimate the term {v, u — Uh). Since dv = dvh and Kv £ 3*^^, we get 
{v, u — Uh) = {Kdvh, d{u — Uh)) — {Kdv, d{u — Uh)) = (dKv, d{u — Uh)) 
= ((/ - TTh)dKv, d{u - Uh)) + {dTThKv, d{u - Uh)) 
= {{I - iih)dKv, d{u - Uh)) + {dP<3*^TThKv, d{u - Uh)). 

Now 

Uh)) < - Uh)\\ < Cr]\\d{u - Uh)\\\\vh\\, 



((/ — 'Kh)dKv, d{u 
while 

{dP<zi-KhKv, d{u - Uh)) = -{d{(T - (Th) 
= -{d{a - ah ) 
< [\\dia-ah)\\ 

But 



{P - Ph)-P<Bi'^hKv - Kv) 
(P - P^hP)^P'Sl'^hKv - Kv) 
-\\p-P^M\]\\P^t^hKv~Kv\, 



\P<ST'!ThKv 



KvW^ = {P^s^nhKv ~ Kv, nhKv - Kv) 

< \\P^-;^nhKv ~ Kv\\\\{I - nh)Kv\\ < \\P<sinhKv - Kv\\5\\v\' 



Hence, 



\\P<s;^TrhKv - Kv\\ <CS\\vh\\. 
Combining these results, and using the previous lemmas, we obtain 

\\vh\\ < Ci7^\\diTrhU-Uh)\\+vMu-Uh)\\+E{u)+S[\\d{a-ah)\\ + \\p-Pf,M\]) 

< C{E{u) + r]E{du) + {rf + 5)[E{da) + E{p)]}. 
The final result of the lemma follows immediately. □ 



It is now an easy matter to prove (39) and so complete the proof Theorem 3.11 
We write 

U-Uh= {u- Phu) + P>8fc(M - Uh) + PjjfcW + P'Sliu - Uh), 



and so (39) follows from Lemmas 3.13 3.14 and 3.16 



To get a feeling for these results, we return to the example mentioned earlier, 
where V,f"^ = Vr+iA''^ ''-{%) and = Vr^^{Th) are used to approximate the 
fc-form Hodge Laplacian with some r > 1. If the domain is convex, we may apply 
elliptic regularity to see that Kf belongs to the Sobolev space H'^h^ for / G L^h^ , 
so dKj E H^\^^^ and d* K f G H^h.^~^ , and then standard approximation theory 
shows that -q = 0{h), S = 0{h^), and fi 
obtain that 



From Theorem 



3.11 



we then 



\\a - ahW + h\\d{a - ah)\\ + h\\u - Uh\\ + h\\p - ph\\ + h'\\d{u - Uh)\\ - Oih'-+') 
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assuming the solution u is sufficiently smooth. That is, all components converge 
with the optimal order possible given the degree of polynomial approximation. 

Finally we note a corollary of Theorem |3.f f | which will be useful in the analysis 
of the eigenvalue problem. 



Corollary 3.17. Under the assumptions of Theorem \3 . 1 1\ there exists a constant 
C such that 

\\K - KhPh\\c(w^.w^) < Cirf + (5 + /i), 
\\dK - dKhPh\\c{w'' ,w>'+^) + \\d*K - d*f^KhPh\\c(w'' ,w''-^) ^ ^V- 
Therefore all three operator norms converge to zero with h. 

Proof. Let / e and set u = Kf, a = d*Kf , p = P?,/, and Uh = KhPhf, 
ah = dlKhPhf- The desired bounds on {K - KhPh)f ^u-ut, d{K - KhPh)f = 



d{u — Uh), and {d*K — d*j^KhPh)f = o" - (ih follow from Theorem 3.11 



since 



E{u)<5\\f\l E{du)+E{a)<r^\\f\l E{da) + E{p) + E{P^u) < \\f\\. 

□ 

We end this section with an estimate of the difference between the true Hodge 
decomposition and the discrete Hodge decomposition of an element of V^^ . While 
this result is not needed in our approach, such an estimate was central to the 
estimation of the eigenvalue error for the Hodge Laplacian using Whitney forms 
made in ^40j . 

Lemma 3.18. Let P'rQU + P^u + PtQ-^u and P<Si^u + Pf,i^u + P'rg'^u denote the Hodge 
and discrete Hodge decompositions of u (£ . Then 

\\P<S'U- P<siu\\ + \\Pf,u-Pi,,u\\ + \\P<su-P<B,u\\ < C{t] + fi)\\u\\v. 

Proof. Let = -Rb*^ and v ~ P'S'P^u'^ — P<s*u. The first estimate follows 
immediately from Lemma |3.12| To obtain the second estimate, we write 

P^U - P^,U = Pf,iP^^U + P^S'U) ~ P^^U = (Pf, - I)P^^U + Pf,{P^'U - P^.u). 

The second estimate now follows directly from ( 28 1 and the first estimate and the 
final estimate follows from the first two by the triangle inequality. □ 

Applied to the case of Whitney forms, both 77 and ^ are 0{h), and so this result 
improves the 0{h\ logh\) estimate of [40, Theorem 2.10]. 

3.6. The eigenvalue problem. The purpose of this section is to study the eigen- 
value problem associated to the abstract Hodge Laplacian ( 18 1. As in the previous 
section, we will assume (W, d) is a Hilbert complex satisfying the compactness prop- 
erty and that the cochain projections tt^ are bounded in £{W'', W''), uniformly in 
h. A pair (A,w) S M x V'', where u 0, is referred to as an eigenvalue/eigenvector 
of the problem (18 1 if there exists (u, p) G V'^^^ x Sj'' such that 

(ct,t) - (dr,u) =0, TeV''-\ 

(42) {da,v) + {du,dv) + {v,p) = X{u,v) veV'', 

In operator terms, u = XKu, a = d*u, p ^ 0. Note that it follows from this system 
that 

x\\ur = \\duf+\\d*ur>o, 
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SO Ku = \^^u. Since the operator K G C{W^ ^W^) is compact and self-adjoint, 
we can conclude that the problem (42 1 has at most a countable set of eigenvalues, 
each of finite multiplicity. We denote these by 

< Ai < A2 < . . . 

where each eigenvalue is repeated according to its multiplicity. Furthermore, when 
W'^ is infinite dimensional, we have lim-,_>oo Aj = 00. We denote by {vi} a corre- 
sponding orthonormal basis of eigenvectors for . 

The corresponding discrete eigenvalue problem takes the form 

K,T)-(dr,u„>-0, r^V;^'\ 

(43) {d(7h,v) + {duh,dv) + {v,ph) = Xh{uh,v), veV/^, 

{uh,q)^o, q<^S)i, 

where Xh G K, and {ah, Uh,Ph) G Vh~^ ^ ^/f ^ -^L "^ith Uh ^ 0. As above, we can 
conclude that p/,, = 0, A/,, > 0, and that A,^^ is an eigenvalue for the operator K^, 
i.e. KhUh — A^^u/j. We denote by 

< Xlji < A2,/!. < . ■ . < ^Nh,h, 




Nh = dimV^', the eigenvalues of problem (43), repeated according to multiplicity, 
and by {wi,^} a corresponding M^'^-orthonormal eigenbasis for V/^. 

Next, we will study how the discrete eigenvalue problem (43 1 converges to the 
eigenvalue problem (42), i.e., how the eigenvalues and eigenvectors of the operator 
Kh converge to those of K. We let Ei and Ei^h denote the 1-dimensional spaces 
spanned by Vi and Viji. For every positive integer j, let m(j) denote the number of 
eigenvalues less than or equal to the jth distinct eigenvalue of the Hodge-Laplace 
problem (42 1. Thus ^^e space spanned by the eigenspaces associated to 

the first j distinct eigenvalues and does not depend on the choice of the eigenbasis. 

The discrete eigenvalue problem (43) is said to converge to the exact eigenvalue 
problem (42) if, for any e > and integer j > 0, there exists a mesh parameter 
ho > such that, for all h < ho, wc have 

(44) max \Xt-Xi,h\<£, and gap 

\i=l 

where the gap between two subspaces E and F of a Hilbert space is defined by 
( 26 1 . The motivation for this rather strict concept of convergence is that it not only 
implies that each eigenvalue is approximated by the appropriate number of discrete 
eigenvalues, counting multiplicities, and that the eigenspace is well-approximated 
by the corresponding discrete eigenspaces, but it also rules out spurious discrete 
eigenvalues and eigenvectors. In particular, this rules out the behavior exemplified 
in Figure |2.6| in Section |2] 

A key result in the perturbation theory of linear operators is that, for eigen- 
value problems of the form we consider, corresponding to the bounded compact 
self-adjoint Hilbert space operators K and Kh, convergence of the eigenvalue ap- 
proximation holds if the operators KhPh converge to K in £{W^, W'^). This result, 
widely used in the theory of mixed finite element eigenvalue approximation, essen- 
tially follows from the contour integral representation of the spectral projection, 
and can be extracted from [BH Chapters HI, IV] or [T3J Section 7]. For a clear 
statement, see [TH]. In fact, as observed in [2D], this operator norm convergence is 
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sufficient as well as necessary for obtaining convergence of the eigenvalue approxi- 
mations in the sense above. As a consequence of Corollary |3.17| we therefore obtain 
the following theorem. 

Theorem 3.19. Let (V, d) be the domain complex of a closed Hilbert complex 
(W,d) satisfying the compactness property, and (V/j,d) a family of subcomplexes 
parametrized by h and admitting uniformly W -bounded cochain projections. Then 
the discrete eigenvalue problems (43l converge to the problem (42l, i.e., (44l holds. 

It is also possible to use the theory developed above to obtain rates of convergence 
for the approximation of eigenvalues, based on the following result (Theorem 7.3 
of |14j). Here, we consider only the case of a simple eigenvalue, but with a small 
modification the results extend to eigenvalues of positive multiplicity. 

Lemma 3.20. If X is a simple eigenvalue and u a normalized eigenvector, then 

\\-^ - < C{\{{K - KhPh)u,u)\ + \\{K - KhPhMl). 

Theorem 3.21. Let (V, d) be the domain complex of a closed Hilbert complex 
{W,d) satisfying the compactness property, and {Vh,d) a family of subcomplexes 
parametrized by h and admitting uniformly W -bounded cochain projections. Let A 
be a simple eigenvalue, u the corresponding eigenvector, and \h be the correspond- 
ing discrete eigenvalue. Let w = Ku, (j) — d*w denote the solution of the mixed 
formulation of the Hodge Laplace problem with source term u, and let Wh = I'ChPhU, 
(j)h — d^Wh denote the corresponding discrete solution. Then 

|A-i - A^^l < C{\\w ~Whfv+U~ + I W - M,w~ Wh)\). 

Proof. To estimate the first term in Lemma |3.20[ we write 

{{K - KhPh)u, u) = {Ku, u) - {KhPnu, P^u) 

— {w, {d*d + dd*)w) — {wh, {d*}^d + dd'^)wh) 

= \\dwf + \\<j>f - \\dWhf - Uhf 

= \\d{w-Wh)\\^ + 2{d{w-Wh),dwh) - U-M^ + 2{(j>-(t)h,<l>) 

= \\diw - Wh)f 2{d{cp - 4>h),WH) + 2(d(0 - 4>h)M 

= \\d{w-WhW - \\<l)-(f>h\? + '^{d{(j)-(l)h),w-Wh), 

where we have used the orthogonality condition 

- 4>h),v) + {d{w - Wh),dv) =0, we Vl: 

in the second last line above. The theorem then follows from this estimate and the 
fact that \\{K-KhPh)u\\l = \\w-Wh\\l. □ 

Order of convergence estimates now follow directly from Theorem |3.11| with / 
replaced by u , u replaced by w, and a replaced by (j). In the example following 
, Vj"-^ = Vr+i^^^^{Th) and Vj^ = VrA''{Th), we saw that 



Theorem 



3.11 



U- 0„|| + h\\d{q^ - + h\\w - w,,\\ + h^\\d{w - w,,)\\ - 0(r+2), 

as long as the domain is convex and the solution w — Ku sufficiently smooth. 
Inserting these results into Theorem 3.21 we find that the eigenvalue error |A— A;i| = 
0{h^'^), double the rate achieved for the source problem. As another example, one 
can check that for the Whitney forms, Vj^'^ = V^K^-^iTh), = V^K^{Th) we 
get |A — Xh\ = 0{h'^), improving on the 0{h\logh\) estimate of [iO] . 
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Remark. It has long been observed that for mixed finite element approximation of 
eigenvalue problems, stability and approximability alone, while sufficient for con- 
vergence of approximations of the source problem, is not sufficient for convergence 
of the eigenvalue problem. An extensive literature developed in order to obtain 
eigenvalue convergence, and a wide variety of additional properties of the finite 
element spaces have been defined and hypothesized. In particular, the first con- 
vergence results for the the important case of electromagnetic eigenvalue problems, 
were obtained by Kikuchi based on the discrete compactness property |65l 119] . and 
the more recent approach by Boffi and collaborators emphasize the Fortid prop- 
erty |18j . In our context, the discrete compactness property says that whenever 
Vh S is uniformly bounded in V'', then there exists a sequence hi converg- 
ing to zero such that u^. converges in W'^, and the Fortid property says that 
lim/j^o 11-^ ~ '^h\\c(v''nv* ,w'') = 0- These additional properties play no role in the 
theory as presented here. We prove eigenvalue convergence under the same sort of 
assumptions we use to establish stability and convergence for the source problem, 
namely subcomplexes, bounded cochain projections, and approximability, but we 
require boundedness in W for eigenvalue convergence, while for stability we only 
require boundedness in V. We believe that subcomplexes with bounded cochain 
projections provide an appropriate framework for the analysis of the eigenvalue 
problem, since, as far as we know, these properties hold in all examples where 
eigenvalue convergence has been obtained by other methods. Moreover, it is easy 
to show that the discrete compactness property and Fortid property hold whenever 
there exist T4^-bounded cochain projections. 

3.6.1. Related eigenvalue problems. Recall that the source problem for the Hodge 



Laplacian, given by (18 1, can be decomposed into the problems (19 1 and (20 1, re- 
ferred to as the *B* and the ?B problem, respectively. More precisely, these problems 
arise if the right hand side / of (18) is in 03^ or S'^, respectively. In a similar way 



the eigenvalue problem (42 1 can be decomposed into a S* and the 05 problem. To 



see this, note if {X,u) g Mx^ is an eigenvalue/eigenvector for (42l, then w G S) 



and therefore u — utQ + uj_, where u<s £ 25 and uj_ G 05^ nV= 3 • It is 
straightforward to check that uj_ satisfies the 05^. eigenvalue problem given by 

(45) (duj_,dv) ^ X{uj_,v), w e S*"^- 

Here uj_ can be equal to zero even if u 7^ is an eigenvector. On the other hand. 



if the pair {X,u±) e M x 3*^^ is an eigenvalue/eigenvector for (45 1, then {X,u±) is 
also an eigenvalue/eigenvector for ( |42| , where a and p are both zero. Furthermore, 
for the extended eigenvalue problem 

(46) {du^,dv) ^ X(uj_,v), veV'', 

where u± is sought in V'^', the nonzero eigenvalues correspond precisely to the 



eigenvalues of ( |45| , while the eigenvalue A = has the eigenspace 3'^- 

The pair {X,u<s) satisfies the corresponding eigenvalue problem given by 

(a,r) -(dT,u<s) =0, tgV''-\ 
(47) . 

{da,v) = X{u<B,v) w e <8 , 

where a = d*u<B € ^l-i ^ V^~'^ = 3^''~^^"''- Furthermore, any solution {X,u<s) 



of (47), where Ufg ^ 0, corresponds to an eigenvalue/eigenvector of the full Hodge 



Laplacian (42). In fact, any eigenvalue of (|42| is an eigenvalue of either the 03* or 
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05 problem, or both, and all eigenvalues of the 58* and 58 problems correspond to 
an eigenvalue of ( |42[). In short, the eigenvalue problem (42 1 can be decomposed into 



the two problems (45 I and (47) 



Also, note that if (A,7if8) is an eigenvalue/eigenvector for the *8'^ problem (47 1 
and a — d*u^ £ 3^*^^^''^, then <t ^ 0, and by taking v = dr in (47l, we obtain 



(48) 



{da.dr) = A((t, r), r e 3 



which is a 58J!_j^ problem. Hence, we conclude that any eigenvalue of the 58 
problem is an eigenvalue of the 58^_j^ problem. The converse is also true. To see 
this, let (A, cr) e M X 3^''"^^^, cr 7^ 0, be a solution of (Esl, and define e 58'' such 
that d*w>B = a. Since is a bijection, this determines u<s uniquely, 

and (AjMig) is an eigenvalue/eigenvector for the 58*"' problem (47 1. Therefore, the 
two problems (47 1 and ( |48[ ) are equivalent, with the same eigenvalues, and with 
eigenvectors related by d*u<;g — a. In particular, the eigenvalues of the full Hodge 
Laplace problem (42) correspond precisely to the eigenvalues of the 58^._]^ and 58^ 
problems. 

In the special case of the de Rham complex, the 58 1 problem is closely related to 
Maxwell's equations. Discretizations of this problem have therefore been intensively 
studied in the literature (e.g., see [12] and the references therein). However, it is 
usually not straightforward to compute the eigenvalues and eigenvectors of the 
discrete 58t problem from the formulation (45) since we do not, a priori, have a 



basis for the corresponding space 3^''" C available. Usually, we will have only 
have a basis for the space Vj^ at our disposal, and a direct computation of a basis for 
3^^ from this basis is costly. A better alternative is therefore to solve the discrete 
version of the extended eigenvalue problem ( 46 1 , and observe that the positive 



eigenvalues, and the corresponding eigenvectors, are precisely the solutions of the 
corresponding problem (45). Alternatively, we can solve the full Hodge Laplace 



problem (43), and then throw away all eigenvalues corresponding to ah 0. 



4. Exterior calculus and the de Rham complex 

We next turn to the most important example of the preceding theory, in which 
the Hilbert complex is the de Rham complex associated to a bounded domain 
in M". We begin by a quick review of basic notions from exterior calculus. Details 
can be found in many references, e.g., [5] [T^ li^ l [5^ [571 [55] . 

4.1. Basic notions from exterior calculus. For a vector space V and a non- 
negative integer k, we denote by Alt'' V the space of real-valued fc-linear forms on 
V. If dimV" 71, then dim Alt '^ V ^ [l). The wedge product A G Alt^+'' V of 
uj e Alt^ V and 77 e Alt'" V is given by 

{u A ?7)(?;i, . . . ,Wj+fc) = ^(signcr)w(w^(i),. . .,i;<,Q))7y(7;^(j+i), . . .,v^(^j+u)), 

a 

where the sum is over all permutations ct of {1, . . . , j + fc}, for which cr(l) < cr(2) < 
• • • < o'(j) and CT(j-l-l) < <T(j + 2) < • • • < (T(j + fc). An inner product on V induces 
an inner product on Alt'' V: 

(^.'7> = 'Y^'^{(^a{\), ■ ■ ■ ■,e-a(k))'n{(^a(\-)-, ■ ■ ■ ,e.a(k)), 7^ € Alt'' V", 
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where the sum is over increasing sequences a : {1, . . . , fc} — s- {1, . . . , n} and ei, . . ., 
e„ is any orthonormal basis (the right-hand side being independent of the choice 
of orthonormal basis). If we orient V by assigning positive orientation to some 
particular ordered basis (thereby assigning a positive or negative orientation to all 
ordered bases, according to the determinant of the change of basis transformation), 
then we may define a unique volume form vol in Alt" V,n = dim V, characterized by 
vol(ei, . . . , e„) = 1 for any positively oriented ordered orthonormal basis Ci, . . . , e„. 
The Hodge star operator is an isometry of Alt'^ V onto Alt""*' V given by 

w A = (★w, /i)vol, u; e Alt'" V, iJ.e Alt""*' V. 

Given a smooth manifold il, possibly with boundary, a differential k-form ui is 
a section of the fc-altcrnating bundle, i.e., a map which assigns to each a; €E an 
element tOx G Alt*' where denotes the tangent space to at x. We write 
C"'A'^(ri) for the space of m times continuously differentiable fc-forms, i.e., forms 
for which 



belongs to C™(fi) for any smooth vector fields Vi. Similarly, we may define C°°A*=, 
(joo^k (gjnooth forms with compact support contained in the interior of fi), etc. If 
Q is a Riemannian manifold, and so has a measure defined on it, we may similarly 
define Lebesgue spaces LpA'' and Sobolev spaces W^A'^ and H"^A'' = Wf'A*^. The 
spaces H™A'^ are Hilbert spaces. In particular, the inner product in L2A*' = H^A'^ 
is given by 



We also write A*^ or A'^{U) for the space of all smooth differential fc-forms, or at 
least sufficiently smooth, as demanded by the context. 

For any smooth manifold Cl and w G A*'(Q), the exterior derivative du) is a 
{k + l)-form, which itself has vanishing exterior derivative: d{duj) = 0. On an 
oriented n-dimensional piecewise smooth manifold, a differential n-form (with, e.g., 
compact support) can be integrated, without recourse to a measure or metric. 

A smooth map (/) : CI' ^ 0. between manifolds, induces a puUback of differential 
forms from to il' . Namely, if w € A'^(ri), the pullback 4)*ijj G A''{Cl') is defined by 

{(p*u})x{vi, ...,Vk) = oj^:^x){D(l)x{vi), . ■ . , D(j)x{vk)) , x G CI' , vi,...,VkG T^^' . 

The pullback respects exterior products and exterior derivatives: 

0* (uJ A T/) = (j)*uj A 4>*ri, (j}*{duj) — dcffuj. 

If (j) is an orientation-preserving diffeomorphism of oriented n-dimensional mani- 
folds, and a; is an n-form on f2, then 



If ri' is a submanifold of fi, then the pullback of the inclusion map is the trace 
map trn n/, written simply tr^/ or tr when the manifolds can be inferred from 
context. We recall the trace theorem which states that if Q.' is a submanifold of 
codimension 1, then the trace map extends to a bounded operator from H^A^{^) to 
L^A'^(f]'), or, more precisely, to a bounded surjection of ff ^A'^(r2) onto ff ^/^A*^(n'). 
A particularly important situation is when Q.' = dCl, in which case Stokes's theorem 



X Wx{Vi{x), . . .,Vn{x)) 





UJ. 
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relates the integrals of the exterior derivative and trace of an (n — l)-form u on an 
oriented n-dimensional manifold-with-boundary fl: 



duj = trw. 

Jn Jao, 



(A common abuse of notation is to write J^^^ lu for the J^^^ tr lu.) Applying Stokes's 
theorem to the differential form oj A rj with uj e A*^^^, 77 e A""*^, and using the 
Leibniz rule d{u} Arj) = dii)Ari+ (— l)'^~^a; A dr], we obtain the integration-by-parts 
formula for differential forms: 



(49) I du) AT] = {-l)^ j u) A d-q+ j trwAtrjj, a; e A'^"^ t? e A" 
in JQ. Jon 



I n—k 

IL LU /\ IL If, LU IV , 1 1 C 1 

Ion 

On an oriented n-dimensional Riemannian manifold, there is a volume form 

vol G A"(fi) which at each point x of the manifold is equal to the volume form on 
T^fl. Consequently, the Hodge star operation takes iv G A'^(r2) to -ku) £ A"~'^(r2) 
satisfying 

for all /Lt € A"~''(r2). Introducing the coderivative operator S : A'' ^ A'^"^ defined 
by 

(50) i^Sco = i-l)''d*uj, 

and setting 77 = the integration- by-parts formula becomes 



(51) {dbj,^i) = {u;,5ijl) + / trwAtr^/x, a; € A" \ /x e A". 

Jaa 

Since the Hodge star operator is smooth and an isometry at every point, every 
property of /c-forms yields a corresponding property for n — k forms. For example, 
the Hodge star operator maps the spaces C""A'', W^A'^, etc., isometrically onto 
(om^n-fe wrn^n-k^ g^^. gy definition, S : A""'' ^ A^-k-i corresponds to d : 
A'^ A'^+^ via the Hodge star isomorphism. Thus for every property of d there is 
a corresponding property of 6. 

In case Q is a domain in M", we may (but usually will not) use the standard 
coordinates of IR" to write a fc-form as 



CO 



tta- dxai A - ■ ■ Adx, 



l<cri<---<crfc<n 



where the € L^(r2) and dxj : K" ^ M is the linear form which associates to a 
vector its jth coordinate. In this case, the exterior derivative is given by the simple 
formula 

d{aa dxai A - ■ ■ A dx^^ ) = ^-^ dxj A dx^i A - ■ ■ A dx^^ , 

extended by linearity to a sum of such terms. If is a domain in M", then a dxi A 
■ ■ -Adxa^), has the value suggested by the notation. In the case of a domain in M", 
the volume form is simply dxi A ■ ■ ■ A dxn- We remark that the exterior derivative 
and integral of differential forms can be computed on arbitrary manifolds from the 
formulas on subdomains on R" and pullbacks through charts. 
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4.2. The de Rham complex as a Hilbert complex. Henceforth we restrict 
attention to the case that is a bounded domain in M" with a piecewise smooth, 
Lipschitz boundary. In this section we show that the de Rham complex is a Hilbert 
complex which satisfies the compactness property, and so the abstract theory of 
Section |3] applies. We then interpret the results in the case of the de Rham complex. 

We begin by indicating how the exterior derivative d can be viewed as a closed 
densely-defined operator from W'' = L'^A''{n) to W^'^^ = L?'K^+^{9?). Let uj G 



Lp'K^iyi). In view of (51 1, we say that rj E L'^ A''^^ (il) is the weak exterior derivative 



of LU if 

The weak exterior derivative of uj, if one exists, is unique and we denote it by duj. 
In analogy with the definition of Sobolev spaces (cf., e.g., [471 Section 5.2.2]), we 
define HA'' to be the space of forms in L^A'^ with a weak derivative in L^A'''^^. 
With the inner product 

this is easily seen to be a Hilbert space and clearly c? is a bounded map from HA'' 
to L^A'^^^. A standard smoothing argument, as in [ITl Section 5.3], implies that 
C°°A'^(n) is dense in HA''. We take HA''{i}) as the domain V'' of the exterior 
derivative, which is thus densely defined in L'^A''{il). Since HA'' is complete, d is 
a closed operator. 

Thus the spaces L^A''{U) and the exterior derivative operators d form a Hilbert 
complex with the associated domain complex 

(52) ^ HA'^ifl) i HA^n) ^ • • • i iJA"(f^) ^ 0. 

This is the de Rham complex. We shall see below that it satisfies the compact- 
ness property. 

To proceed, we need to identify the adjoint operator d* and its domain . 
Using the surjectivity of the trace operator from ff^A'(ri) onto H^^^ A'{d^) and 



the integration- by-parts formula (49 1, we can show that (cf., [8, page 19]) the trace 



operator extends boundedly from HA''{^) to H ^/^A''(9ri), and that (49 1 holds 



for u) G HA'' ^, 1] G H^A^ '' where the integral on the boundary is interpreted 
via the pairing of H-^/'^{dn) and H^^^{dn). Equivalently, we have an extended 



version of (51 ) 



(53) {dLu,^j.) = {uj,S^i) + trtjAtr*/i, lu G HA''~'^ , e H^ A'' . 

Jan 

Of course, there is a corresponding result obtained by the Hodge star isomorphism 



which interchanges d and S. After reindexing, this is nothing but the fact that (53) 
holds also for a; e H^A''^^, p e H*A'' where 

(54) iJ*A'= :=*(iJA"-'^). 

Note that H* A'' consists of those differential form in L'^A'' for which a weak 
coderivative exists in L'^A''^^, where the weak exterior coderivative is defined in 
exact analogy to the weak exterior derivative. Its inner product is 

The space H*A'' is isometric to HA"^'' via the Hodge star, but is quite different 
from HA''. 
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We also make use of the trace defined on HK'' to define the subspace with 
vanishing trace: 

HK^{n) = {uj£ HA'^in) I traow = 0}. 
Correspondingly, for cj G H* , the quantity trT^cj is well defined, and we have 

(55) H*h^{Q.) ★iJA""'^ = G H*K^{Q.) \ trgn = 0}. 



From ( 53 1 , we have 



(56) {duj,^i) ^ {lj,5^i), uj e HK^-^, ^le H*K^. 

(We certainly have (56 1 with the stronger condition lo E H^A'^^^, but then we can 



extend to all uj £ HAF ^ by continuity and density.) Of course, the corresponding 
result, where uj e HK^'^, ^ e H*A'', holds as well. 

Theorem 4.1. Let d be the exterior derivative viewed as an unbounded operator 
^2^A;-i _^ 7j2^fe domain HA''. Then the adjoint d* , as an unbounded operator 
L^A^ L'^A^^^, has as its domain H*A^, and coincides with the operator 6 defined 



m (501 



Proof. We must show that for /i G L^A'^, there exists oj G Li^A^ ^ such that 
(57) {^jL.dv) ^ veHA^^^, 

if and only if /i e H*A^ and uj = Sji. The if direction is immediate from (56 1. 
Conversely, if (57 1 holds, then /i has a weak exterior coderivative in , namely 



5/i = u). Thus /i e H* A . Integrating by parts we have 

tri/ Atr*^ = {^i,dv) - (uj,iy) =0, v E iJ^A''"^ 



Ian 

which implies that tr*fi — 0, i.e., /i G H* A^ . □ 

As a corollary, we obtain a concrete characterization of the harmonic forms: 

(58) S)^ = {lo E HA^ C]H*A^\duj ^Q,5uj = Q}. 

In other words, a fc-form is harmonic if it satisfies the differential equations duj = Q 
and 5uj = Q together with the boundary conditions trik-w = 0. 

If the boundary of il is smooth, then HA^ n H*A*' is contained in H^A^ [5^ . 
and hence, by the Rellich theorem, we obtain the compactness property of Sec- 
tion 3.1.3| For a general Lipschitz boundary, e.g., for a polygonal domain, the 
inclusion of HA'^ n H*A^ in H^A'' need not hold, but the compactness property 
remains valid [75]. Thus all the results of Section |3] apply to the de Rham complex. 
In particular, we have the Hodge decomposition of L^A'^ and of HA'', the Poincare 
inequality, well-posedness of the mixed formulation of the Hodge Laplacian, and all 
the approximation results established in Sections |3. 3^3.61 We now interpret these 
results a bit more concretely in the present setting. 



First of all, the cohomology groups associated to the complex ( 52 1 are the 
de Rham cohomology groups, whose dimensions are the Betti numbers of the do- 
main. Turning next to the Hodge Laplacian problem given in the abstract case by 



( [T8| , we get that {a, u,p) G HA''^^ x HA'^ x is a solution if and only if 

(59) a ~ Su, da + 6du = f — p in 51, 

(60) tr-ku — 0, trirdu — on dfl, 

(61) U-LiD^ 
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The first differential equation and the first boundary condition are implied by the 



first equation in (18 1, and the second differential equation and boundary condi- 



tion by the second, while the third equation in (18 1 is simply the side condition 
u _L i]*^. Note that both boundary conditions are natural in this variational for- 
mulation: they are implied but not imposed in the spaces where the solution is 
sought. Essential boundary conditions could be imposed instead. We discuss this 
in Section [62] 

To make things more concrete, we now restrict to a domain 17 C M"^, and consider 
the Hodge Laplacian for fc- forms, k — 0,1,2, and 3. We also discuss the 23* and ?B 



problems given by (19 1 and (20 1 for each k. We shall encounter many of the most 
important partial differential equations of mathematical physics: the Laplacian, the 
vector Laplacian, div-curl problems, and curl-curl problems. These PDEs arise in 
manifold applications in electromagnetism, solid mechanics, fluid mechanics, and 
many other fields. 

We begin by noting that, on any oriented Riemannian manifold of dimension n, 
we have a natural way to view 0-forms and n-forms as real-valued functions, and 
1-forms and (n— l)-forms as vector fields. In fact, 0-forms are real- valued functions 
and 1-forms are covector fields, which can be identified with vector fields via the 
inner product. The Hodge star operation then carries these identifications to n- 
forms and {n — l)-forms. In the case of a 3-dimensional domain in M"^, via these 
identifications all fc-forms can be viewed as either scalar or vector fields (sometimes 
called proxy fields). With these identifications, the Hodge star operation becomes 
trivial in the sense if a certain vector field is the proxy for, e.g., a 1-form lo, the 
exact same vector field is the proxy for the 2-form -kuj. Via proxy fields, the exterior 
derivatives coincide with standard differential operators of calculus: 

d° ~ grad, ~ curl, d^ — div, 

and the de Rham complex ( [52| is realized as 

^ H\n) ^ if (curl; n) ^ H{dW; n) ^ L^{n) ^ 0, 

where 

if (curl; n) ^ {u:n^R^\ue L^, curlu e }, 
H{div; = {u:fl^R^\ue L^,dWue L^}. 

The exterior coderivatives S become, of course, — div, curl, and — grad, when acting 
on 1-forms, 2-forms, and 3-forms, respectively. The trace operation on 0-forms is 
just the restriction to the boundary, and the trace operator on 3-forms vanishes 
(since there are no nonzero 3-forms on The trace operator from 1-forms on 

f2 to 1-forms on the boundary takes a vector field u on f2 to a tangential vector 
field on the boundary, namely at each boundary point x, (trM)^, is the tangential 
projection of u^- For a 2-form u, the trace corresponds to the scalar u ■ n (with n 
the unit normal) at each boundary point. 

4.2.1. The Hodge Laplacian for k ~ 0. For fc = 0, the boundary value problem 



(59l-(61) is the Neumann problem for the ordinary scalar Laplacian. The space 
HA^^ is understood to be 0, so a vanishes. The harmonic form space Sj'^ consists 
of the constant functions (we assume is connected — otherwise io" would consist 
of functions which are constant on each connected component), and p is just the 
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average of /. The first differential equation of (59 1 vanishes, and the second gives 
Poisson's equation 

— div grad u = f ~ p in fi. 



Similarly, the first boundary condition in (60 1 vanishes, while the second is the 
Neumann condition 

gradw • n = on dCl. 



The side condition (61 1 specifies a unique solution by requiring its average value to 
be zero. 

Nothing additional is obtained by considering the split into the 05* and S sub- 



problems, since the latter is trivial. Furthermore, the eigenvalue problem (42 1 is 
precisely the corresponding eigenvalue problem for the scalar Laplacian, with the 
eigenspace filtered out. 

4.2.2. The Hodge Laplacian for k — 1. In this case the differential equations and 
boundary conditions are 

(T = — div u, grad a + curl curl u — f — p in fi, 
(62) J f ^ 

u ■ n = 0, curlw x n — on dil, 

which is a formulation of a boundary value problem for the vector Laplacian 
curl curl — grad div. (Here we have used the fact that the vanishing of the tan- 
gential component of a vector is equivalent to the vanishing of its cross product 
with the normal.) The solution is determined uniquely by the additional condi- 
tion that u be orthogonal to Sj^ , which in this case consists of those vector fields 
satisfying 

curlp = 0, div p =0 in fi, p ■ n = 0, on dfl. 

The dimension of Sj-^ is equal to the first Betti number, i.e., the number of handles, 
of the domain (so — if the domain is simply-connected). 



The OS* problem (19) is defined for vector fields / which are orthogonal to 



both gradients and the vector fields in Sj^. In that case, the solution to (62 1 has 
cT = and p = 0, while u satisfies 

curl curl u = f, div u = Q infi, u ■ n — Q, curl u x ri = on dU,. 

The orthogonality condition u _L again determines the solution uniquely. 

Next we turn to the problem. For source functions of the form / — grad F for 



some F g , which we normalize so that J^F = 0, (62 1 reduces to the problem 
of finding a G and u S 53^ = gradH^ such that: 

(T = — div u, grad cr = f in 51, u • n = 0, on dfl. 

The differential equations may be simplified to — grad div u = /, and the condition 
that u G 33^ can be replaced by the differential equation curlu = 0, together with 
orthogonality to Sj^. Now grad((T — F) = and J^cr = — /^^ u ■ n = = JqF, 
a = F, and we may rewrite the system as 

— div u — F, curl u = in 17, u • n = 0, on dfl, 

which, again, has a unique solution subject to orthogonality to Sj^. 



so 



The eigenvalue problem ( 42 1 is the corresponding eigenvalue problem for the 



vector Laplacian with the boundary conditions of ( 62 1 , and with the eigenspace 



Sj^ of the eigenvalue A = filtered out. As mentioned above, in this case the 
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eigenvalue problem, given by (451, is important for models based on Maxwell's 
equations. This problem takes the form 

(63) curlcurlu — Xu, divu — in f2, u ■ n — 0, curlw x rt = on dfl, u _L Sj^. 

4.2.3. The Hodge Laplacian for k — 2. The differential equations and boundary 
conditions are 

a = curl u, curl a — grad div u = f — p in fi, 
^^^^ u X n — 0, divu — on dCl. 

This is again a formulation of a boundary value problem for the vector Laplacian 
curl curl — grad div, but with different boundary conditions than for (62 1, and this 
time stated in terms of two vector variables, rather than one vector and one scalar. 
This time uniqueness is obtained by imposing orthogonality to Sj^, the space of 
vector fields satisfying 

curlp = 0, divp = in fi, p x n = 0, on dfl, 

which has dimension equal to the second Betti number, i.e., the number of voids in 
the domain. 

The problem arises for source functions of the form / = grad F for some 
F e {H^ is the function space corresponding to H*A^). We find (7 = 0, and u 
solves 

— div u ~ F, curlu = in fi, u x rt = 0, on dCl, 

i.e., the same differential equation as for *B^, but with different boundary condi- 
tions, an extra assumption on F, and, of course, now uniqueness is determined by 
orthogonality to 

If div/ — and / _L i^^, we get the 03^ problem for which the differential 
equations are a = curlu, curler — f and the condition divu = arising from the 
membership of u in *B^. Thus u solves 

curl curl u = /, divu = in fi, u x n = 0, on dfl, 

the same differential equation as for 03^, but with different boundary conditions. 

The eigenvalue problem ( 42 1 is the corresponding eigenvalue problem for the 
vector Laplacian with the boundary conditions of (64), and with the eigenspace Sj^ 
of the eigenvalue A = filtered out, while the corresponding 05^ problem, of the 
form ( 47 1 , takes the form 

curl curl ?i = Au, divu = in fl, u x n = 0, on dfl, u _L S^^ . 



Note that this is the same problem as the OS^^ eigenvalue problem(63l, but with 
different boundary conditions. However, if we define a = curlu, then it is straight- 
forward to check that the pair (A,(t) will indeed solve the ^\ problem (63 1. This 
is an instance of the general equivalence between the S*^ problem, given by (47 1, 
and the corresponding 03^. -^ problem, which was pointed out in Section 



3.6.1 



4.2.4. The Hodge Laplacian for fc = 3. In this case the Hodge Laplacian problem, 
which coincides with the 05'^ problem, is 

(J — — grad u, div a — f in U,, u = on 9f2, 

which is the Dirichlet problem for Poisson's equation. There are no non-zero har- 
monic forms, and the problem has a unique solution. Furthermore, the eigenvalue 
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problem (42 1 is the corresponding eigenvalue problem for the scalar Laplacian with 



Dirichlet boundary conditions. 

5. Finite element approximation of the de Rham complex 

Our goal in this section is to discretize the de Rham complex so that we may 
apply the abstract results on approximation of Hilbert complexes from Section [3] 
Hence we need to construct finite dimensional subspaces of HA''{U). As we 
saw in Section |3.4| the key properties these spaces must possess is, first, that 
dA^ C A^"'"^ so they form a subcomplex (Ah, d) of the de Rham complex, second, 
that there exist uniformly bounded cochain projections tt/j from (L^A, d) to (A?i, d), 
and third, good approximation properties. We may then use these spaces in a 



Galerkin method based on the mixed formulation, as described in Section 3.4 
and the error estimates given in Theorems |3.9| and |3.11| bound the error in the 
Galerkin solution in terms of the approximation error afforded by the subspaces. In 
this section we will construct the spaces A^ as spaces of finite element differential 
forms and show that they satisfy all these requirements, and can be efficiently 
implemented. 

As discussed in Section [2] a finite element space is a space of piecewise poly- 
nomials which is specified by the finite element assembly process, i.e., by giving 
a triangulation of the domain, a finite dimensional space of polynomial functions 
(or, in our case, differential forms) on each element of the triangulation, called the 
shape functions, and a set of degrees of freedom for the shape functions associated 
with the faces of various dimensions of the elements, which will be used to deter- 
mine the degree of interelement continuity. Spaces constructed in this way can be 
implemented very efficiently, since they admit basis functions with small support 
and so lead to sparse algebraic systems. 

In Section |5.1| we will discuss the spaces of polynomial differential forms which 
we will use as shape functions. In Section |5.2[ we specify degrees of freedom for 
these spaces on simplices, and study the resulting finite element spaces of differential 
forms. In particular, we show that they can be combined in a variety of ways to 
form subcomplexes of the de Rham complex. In Section |5.3| we briefly describe 
recent work related to the implementation of such finite elements with explicit local 



bases. In Sections 5.4 and 5.5 we construct L^-bounded cochain projections into 
these spaces and obtain error estimates for them. Finally Section |5.6| is simply a 
matter of collecting our results in order to obtain error estimates for the resulting 
approximations of the Hodge Laplacian. Many of the results of this section have 
appeared previously, primarily in and therefore many proofs are omitted. 

5.1. Polynomial differential forms and the Koszul complex. In this section, 
we consider spaces of polynomial differential forms, which lead to a variety of sub- 
complexes of the de Rham complex. These will be used in the next section to con- 
struct finite element spaces of differential forms. The simplest spaces of polynomial 
differential fc- forms are the spaces 7'^A''(R") consisting of all differential fc-forms 
on M" whose coefficients are polynomials of degree at most r. In addition to these 
spaces, we will use another family of polynomial form spaces, denoted 7^~A'^'(M"), 
which will be constructed and analyzed using the Koszul differential and the asso- 
ciated Koszul complex. Spaces taken from these two families can be combined into 
polynomial subcomplexes of the de Rham complex in numerous ways (there are 
essentially 2"^^ such subcomplexes associated to each polynomial degree). These 
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will lead to finite element de Rham subcomplexes, presented in Section |5. 2 [ Some 
of these have appeared in the literature previously, with the systematic derivation 
of all of them first appearing in [6]. 

5.1.1. Polynomial differential forms. Let and 7i,.(IR") denote the spaces of 
polynomials in n variables of degree at most r and of homogeneous polynomial 
functions of degree r, respectively. We interpret these spaces to be the zero space 
if r < 0. We can then define spaces of polynomial differential forms, 
7ij.A'^(]R"), etc., as those differential forms which, when applied to a constant vector 
field, have the indicated polynomial dependence. For brevity, we will at times 
suppress E" from the notation and write simply Vr, "Hr, Vr^^ , etc. 

The dimensions of these spaces are easily calculated: 

dimP^A'^XlR") = dimT'rlM") • dim Alt M" 

/ n + r\ / n\ / r + k\ / n + r 
\ n J \k J \ r J \n — k 

and dimH^A'=(M") = dimP^A^^XM""!). 

For each polynomial degree r > 0, we get a homogeneous polynomial subcomplex 
of the de Rham complex: 

(65) ^ HrA" ^ Hr^iA^ ^ . . . ^ Hr-nA'' -> 0. 

We shall verify below the exactness of this sequence. More precisely, the coho- 
mology vanishes if r > and also for r = except in the lowest degree, where 
the cohomology space is M (refiecting the fact that the constants are killed by the 
gradient). 

Taking the direct sum of the homogeneous polynomial de Rham complexes over 
all polynomial degrees gives the polynomial de Rham complex: 

(66) ^ VrA" i Vr-lA^ i . . . ^ Vr-nA'' -> 0, 

for which the cohomology space is M in the lowest degree, and vanishes otherwise. 

5.1.2. The Koszul complex. Let x £ E". Since there is a natural identification of E" 
with the tangent space TqE" at the origin, there is a vector in TqE" corresponding 
to X. (The origin is chosen for convenience here, but we could use any other point 
instead.) Then the translation map ?/ 1— > y + a; induces an isomorphism from T^K^ 
to Tj-E", and so there is an element X{x) e T^W^ corresponding to x. (Essentially 
X{x) is the vector based at x which points opposite to the origin and whose length is 
|a:|.) Contraction with the vector field X defines a map k, from A''"(E") to A'''^^(E") 
called the Koszul differential: 

{kUj)x{vi, . . .,Vk^l) = UJx[X{x),Vl,.. . ,Wfc-l). 

It is easy to see that k is a graded differential, i.e., 

KO K = 

and 

k{u; a 77) = (kuj) a r/ + (-l)''^ A (nrj), co £ A'', rj £ A'. 
In terms of coordinates, if tjj. = a{x) dx^-^ A • • • A dx„^ , then 

k 

{k.uj)x = ^(-l)*+^a(x) x„(i) dx„^ A • • • A dx„(i) A • • • A dx^^, 

i=l 



50 



DOUGLAS N. ARNOLD, RICHARD S. FALK, AND RAGNAR WINTHER 



where the notation dxa^ means that the term is omitted. Note that k maps Ti^A*^ 
to Tir+iA*^"^, i.e., k increases polynomial degree and decreases form degree, the 
exact opposite of the exterior derivative d. 

The Koszul differential gives rise to the homogeneous Koszul complex ^Mi Chap- 
ter 3.4.6], 

(67) ^ 7i,_„A" ^ 7i,_„+iA"-i ^ . . . ^ Hrh^ ^ 0. 

We show below that this complex is exact for r > 0. Adding over polynomial 
degrees, we obtain the Koszul complex (for any r > 0), 

^ n.-„A" ^ Vr-n+l^''-^ ^ . . . ^ T'.AO ^ 0, 

for which all the cohomology spaces vanish, except the rightmost, which is equal to 
M. 

To prove the exactness of the homogeneous polynomial de Rham and Koszul 
complexes, we note a key connection between the exterior derivative and the Koszul 
differential. In the language of homological algebra, this says that the Koszul 
differential is a contracting homotopy for the homogeneous polynomial de Rham 
complex. 

Theorem 5.1. 

(68) {dn + Kd)uj = {r + k)ijj , w G TY^A''. 

Proof. This result can be established by a direct computation or by using the 
homotopy formula of differential geometry (Cartan's magic formula). See [8^ for 
the details. □ 



As a simple consequence of ( 68 1 , we prove the injectivity of d on the range of k 
and vice versa. 

Theorem 5.2. // dnuj — for some u G VA, then ku) = 0. // ndw = for some 
to G VA, then duj — 0. 

Proof. We may assume that uj G JirA^ for some r, fc > 0. If r = fc = 0, the result 
is trivial, so we may assume that r + k > Q. Then (r + k)mjj = n^dn + Kd)uj = 0, 
if dnu) — 0, so KUJ = in this case. Similarly, (r + k)doj = d{dK + Kd)uj — 0, if 
Hiduj = 0. □ 



Another easy application of ( 68 1 is to establish the claimed cohomology of the 
Koszul complex and polynomial de Rham complex. Suppose that u) G TirA^ for 



some r. A: > with r + fc > 0, and that kuj = 0. From ( 68 ) , we see that ui = kt] with 



rj — du)/{r + k) G Tir-iA^^^ . This establishes the exactness of the homogeneous 



Koszul complex ( 67 1 (except when r = and the sequence reduces to ^ M 0) . 



A similar argument establishes the exactness of ( 65 1 . 

Another immediate but important consequence of ( 68 1 is a direct sum decom- 
position of TirA^ for r, fc > with r + fc > 0: 

(69) UrA'' = tiHr-iA^"^^ e d'Hr+iA''-'^ . 

Indeed, if w G HrA^, then rj = duj/{r + fc) G 7Yr-iA''+^ and /i = Kuj/{r -I- fc) G 
Hr+iA''-^ and cj = kt? + d^x, so HrA^ = kH^-iA'^+i + dUr+iA''-'^. Also, if 
Lo G KHr-iA^^^ n dTir+iA^^^ , then duj — klu = Q (since dod=KOK = Q), and so. 
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by (|68]), w = 0. This shows that the sum is direct. Since Vr^^ = ©j^o^j^''' 
also have 

The exactness of the Koszul complex can be used to compute the dimension of 
the summands in (69) (c.f. [8J). 

Theorem 5.3. Let r > 0, I < k < n, for integers r, k, and n. Then 



(70) dimKH,.A'=(M") = dimd7i:^+iA'^-i(M") = 



n + r\ /r + fc— 1 



- k 



k-l 



5.1.3. The space V^A''. Let r > 1. Obviously, VrA'' = Vr-iA^ + 'HrA^. In view 
of (69 1, we may define a space of fc-forms intermediate between Vr-iA'^ and VrA^ 

by 

V-A'' = Vr-lA^ + KHr-lA''+^ = Vr-lA'' + nVr-lA'^+K 

Note that the first sum is direct, while the second need not be. An equivalent 
definition is 

V-A'' = {uje VrA^ I Kuj e VrA^-^ }. 

Note that p- A^ = VrA° and T'-A" = Pr-iA", but for < fc < n, V~A'' is 
contained strictly between Vj.-iA'' and VrA^. For r < 0, we set V^A'' — 0. 
From ( 70 1 , we have 

dimP,7A''(M") = dimVr-iA^ + dim nHr-iA^^^ 



n + r — 1 
n 

r + fc - 1 
k 

where the last step is a simple identity. 



n 
k 

n + r ' 
n ~ k I ' 



n + r — 1 
n — fc — 1 



r + fc - 1 
fc 



We also note the following simple consequences of Lemma 5.2 



Theorem 5.4. If u) e V^. A^ and duj = 0, then uj e T'r-iA . Moreover, dVrA'^ = 
dV-A^. 

Proof. Write uj = uji + kcj2 with ui G Vr-iA^ and uj2 G 7'r-iA'''+^ Then 

= dKUJ2 = => KLJ2 = UJ e Vr-lA'', 

showing the first result. For the second it suffices to note that VrA'^ = A^ + 
dP^+iA*-!. □ 



Remark. We defined the Koszul differential as contraction with the vector field X, 
where X{x) is the translation to x of the vector pointing from the origin in M" to 
X. The choice of the origin as a base point is arbitrary — any point in M" could be 
used. That is, if y G M", we can define a vector field Xy by assigning to each point 
X the translation to x of the vector pointing from y to x, and then define a Koszul 
differential Ky by contraction with Xy. It is easy to check that for oj £ Vr-iA'^'^^ 
and any two points y,y' E M", the difference KyUj — Kyru e T'^-iA'^. Hence the 
space 

V-A'' = Vr-lA'' + KyVr-lA''+^ 
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does not depend on the particular choice of the point y. This observation is impor- 
tant, because it allows us to define 'P^A'^{V) for any affine subspace V of M". We 
simply set 

p-A'=(v) - n-iA'=(y) + KyVr-iA''+\v), 

where y is any point of V. Note that if a; G 7'^A'''(M"), then the trace of a; on F 
belongs to V-A'^iV). 

Remark. The spaces 7'^A'^(IR") are affine-invariant, i.e., if : M" — > M" is an afRne 
map, then the puUback 0* maps this space into itself. Of course, the full polyno- 
mial space VrA^{M.") is affine-invariant as well. In [ff, Section 3.4], all the finite 
dimensional affine-invariant spaces of polynomial differential forms are determined. 
These are precisely the spaces in the V and families together with one further 
family of spaces which is of less interest to us. 



5.1.4. Exact sequences of polynomial differential forms. We have seen the polyno- 



mial de Rham complex ( 66 1 is a subcomplex of the de Rham complex on M" for 
which cohomology vanishes except for the constants at the lowest order. In other 
words, for any r > the sequence 

(71) M ^ VrA^ ^ Vr-lA^ ^ • • • ^ Vr-nA"" ^ 

is exact, i.e., is a resolution of M. 

As we shall soon verify, the complex 

(72) M ^ V-A° ^ V~A^ ^ ...^ p-A" ^ 

is another resolution of M, for any r > 0. Note that in this complex, involving the 



A spaces, the polynomial degree r is held fixed, while in (71 1, the polynomial 



degree decreases as the form order increases. Recall that the 0th order spaces 



in these complexes, VrA^ and V^. A^ , coincide. In fact, the complex ([tT]) is a 



subcomplex of (72), and these two are the extreme cases of a set of 2'^~^ different 
resolutions of M, each a subcomplex of the next, and all of which have the space VrA^ 
in the 0th order. To exhibit them, we begin with the inclusion K ^ VrA^ . We may 
continue the complex with either the map d : VrA^ Vr-iA^ or d : VrA^ V~A^, 
the former being a subcomplex of the latter. With either choice, the cohomology 
vanishes at the first position. Next, if we made the first choice, we can continue 
the complex with either d : Vr-iA^ 'Pr-2A'^ or d : Vr-iA^ T^r-i^^- 
if we made the second choice, we can continue with either d : A^ T'r-iA^ 
or d : V^A^ V^A^. In the first case, we may use the exactness of (71) to 



see that the second cohomology space vanishes. In the second case, this follows 



from Lemma 5.4 Continuing in this way at each order, /c = 1, . . . , n — 1, we have 
two choices for the space of fc-forms (but only one choice for k = n, since T'j.-iA" 
coincides with 'P~A"), and so we obtain 2"~^ complexes. These form a totally 
ordered set with respect to subcomplex inclusion. For r > n these are all distinct 
(but for small r some coincide because the later spaces vanish). 
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In the case n = 3, the four complexes so obtained are: 

M ^ PrA° ^ Pr-lA^ ^ Vr-2A'^ 

R ^ Vrh° -^-^ Vr-l^^ ^r"-!^' 

M ^ VrA^ — ^ V-A^ — ^ Vr-lA^ 
R^VrA^ — ^ V-A^ — ^ p-A2 

5.2. Degrees of freedom and finite element differential forms. Having in- 
troduced the spaces of polynomial differential forms 7',.A'^(M") and 7'^A'^(M"), we 
now wish to create finite clement spaces of differential forms. Wc begin with nota- 
tion for spaces of polynomial differential forms on simplices. If / is a simplex (of 
any dimension) in R", we define 

VrA'^if) = trRnjP,A'=(M"), VrA^if) = {u e VrA^U) I tr/,o/a; = 0}. 

The spaces P^A'^(/) and V~A''{f) are defined similarly. 

Now let Q be a bounded polyhedral domain which is triangulated, i.e., parti- 
tioned into a finite set T of n-simpliccs determining a simplicial decomposition 
of ri. This means that the union of the elements of T is the closure of fi, and 
the intersection of any two is either empty or a common subsimplex of each. By 
way of notation, for any simplex T we denote by Ad{T) the set of subsimplices 
of T of dimension d, and by A(T) the set of all subsimplices of T. We also set 
Ad(T) := {jrerMT) and A(T) := Utst ^(2^)- 

Corresponding to the V and V~ families of spaces of polynomial differential 
forms, we will define two families of spaces of finite element differential forms with 
respect to the triangulation T, denoted VrA''{T) and V~A''{T). We shall show 
that these arc subspaccs of HA'^{il) and can be collected, in various ways, into 
subcomplexes of the de Rham complex. 

The spaces VrA'^{T) and ■p~A'^(T) will be obtained by the finite element as- 
sembly process. For each T G T, wc will choose the corresponding polynomial 
space VrA'^iT) or 'P~A'^{T) to be used as shape functions. The other ingredient 
needed to define the finite element space is a set of degrees of freedom for the shape 
function spaces, that is, a basis for the dual space, in which each degree of freedom 
is associated with a particular subsimplex. When a subsimplex is shared by more 
than one simplex in the triangulation, we will insist that the degrees of freedom 
associated with that subsimplex be single-valued in a sense made precise below, 
and this will determine the interelement continuity. 

The degrees of freedom for VrA'^iT) and V~A'^{T) which we shall associate to 
a d-dimcnsional subsimplex / of T will be of the following form: for some (d — k)- 
form r] on /, the functional will be a; i-^- try,/ oj Arj. The span of all the degrees 
of freedom associated to / is a subspace of the dual space of the shape function 
space, and so we obtain a decomposition of the dual space of the shape functions 
on T into a direct sum of subspaces indexed by the subsimplices of T. It is really 
this geometric decomposition of the dual space that determines the interelement 
continuity rather than the particular choice of degrees of freedom, since wc may 
choose any convenient basis for each space in the decomposition and obtain the 
same assembled finite element space. 







^0, 


Vr- 


.2A3- 




Vr- 


.2A3- 




Vr- 


-iA3- 


0. 
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The geometric decompositions of the dual spaces of Vr^^{T) or V,. K^{T) are 
given specifically in the following theorem, which is proven in [8, Sections 4.5 and 
4.6]. 

Theorem 5.5. Let r, k, and n be integers with < k < n and r > 0, and let T be 

an n-simplex in M". 

1. To each f G A(r), associate a space W^{T, /) C Vr^^{T)*: 



Then W^{T,f) — 7',^j._(jijjj^A'^'™^"'^(/) via the obvious correspondence, and 

Vrh^\TY^ W^{T,J). 
/eA(T) 

2. To each f e A(r), associate a space W^-{T, f) C 'P~A''{T)*: 

(T, /) = I CO trrj Vr+k-dnn /-i A^^™ -^^^ (/) | . 

Then W^^(T, /) = 7-'.r+fc-dim/-iA'^'"^-^^'^(/) via the obvious correspondence, and 

v-K\Ty^ wi:-{T,f). 

feA(T) 

Note that the spaces W^{T,f) and W}^~{T,f) vanish if dim/ < k. Note also 
that the dual space of Vr^^{T) is expressed in terms of spaces in the family, and 
vice versa. This intimate connection between the V and families of spaces of 
polynomial differential forms is most clearly seen in the following algebraic propo- 
sition, which is closely related to Theorem 5.5 and is also proved in Sections 4.5 
and 4.6]. 

Lemma 5.6. With r, k, n as above, 

Vr^\Ty - V;^^^^^K''-\T) and -p-A^T)* - Vr+k-n-iA''~HT). 



With the decompositions in Theorem |5.5| we can define the finite element spaces. 
Thus VrA^iT) consists of all forms w G L?K^{^) such that wIt belongs to the shape 
function space VrA^{T) for all T T, and for which the quantities try w A 77 are 
single-valued for all / £ A(T) and all 77 G 'P~^f._^^^j:A'^™^~''{f). More precisely, 
this means that if / is a common face of Ti,T2 G T, then 

trTi,/(co|Ti) = trT2,/(w|Tj A 7] 

for all such / and 77. The family of spaces is defined analogously. 

The degrees of freedom determine the amount of interelement continuity enforced 
on the finite element space. Of course we need to know that the assembled spaces 
belong to HK^{Q). In fact, the degrees of freedom we imposed enforce exactly 
the continuity needed, as shown in the following theorem. This is proved in [8j 
Section 5.1], where the equations below are taken as definitions, and it is shown 
that the assembly process leads to the same spaces. 
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Theorem 5.7. 

VA^{T) = {w e i/A'=(r!) I lo\t e P,A^ T g T}, 
V^K'^iT) = {w e i?A'=(r!) I cj|t e A^ T e T}. 
Next we note that these spaces of finite element difi^erential forms can be col- 



we 



lected in subcomplexes of the de Rham complex. In view of Theorem 5.7 
have dTrA'^iT) C Vr-iA''+'^{T) and dV-K ^{T) C V-K^+^{T). Corresponding 



to the resolutions defined in Section 5.1.4[ we obtain 2" ^ de Rham subcom- 



plexes for each value of r and each mesh T. Each complex begins with the space 
VrA'^{T) = 'P~hP{T). The maps making up the subcomplexes are all one of the 
following types: 

(73) \ or I — ^ \ or 

{v;+,K^-\r)\ { VsAHV 

for some s and some fc. In Section |5.5| we shall show that all the subcomplexes 
admit bounded cochain projections. These 2"^^ complexes, which are distinct for 
r > n, are linearly ordered by inclusion. The maximal complex is 

^ V-A°iT) ^ r-A\T) ^ • • • ^ V^K'^iT) ^ 0. 

The spaces K^{T) in this complex are referred to as the higher order Whitney 
forms, since for r = 1, this is exactly the complex introduced by Whitney [88 . The 
minimal complex with this starting space is the complex of polynomial differential 
forms 

^ VrA°{T) ^ Vr-iA\T) ^ • • • i 7'^_„A"(T) ^ 0. 

This complex was used extensively by Sullivan [83l [84] and is sometimes referred 
to as the complex of Sullivan- Whitney forms |15| . It was introduced into the finite 
element literature in |37j . The intermediate complexes involve both higher order 
Whitney spaces and full polynomial spaces. 

Finally, we note that the degrees of freedom used to define the space determine 
a canonical projection 1^ : CA''{Q) VrA^iT). Namely, Ih^j G VrA^iT) is 
determined by 

trrjiu^ - 4c.) A 77 = 0, t? G 7'-+,_,,„^Ad'--^-'=(/), / G A(r). 

Similar considerations apply to 'P~ A^ {T) . The canonical projection can be viewed 
as a map from the (sufficiently) smooth de Rham complex to any one of the com- 



plexes built from the maps ( 73 1 . The degrees of freedom were chosen exactly so that 
the canonical projection is a cochain map, i.e., commutes with d. This is verified 
using Stokes's theorem. See [H Theorem 5.2]. 

Remark. Given a simplicial triangulation T, we have defined, for every form degree 
k and every polynomial degree r, two finite element spaces of differential fc- forms 
of degree at most r: VrA^{T) and V~ A''{T). In n = 2 and n = 3 dimensions, 
we may use proxy fields to identify these spaces of finite element differential forms 
with finite element spaces of scalar and vector functions. The space VrA^iT) cor- 
responds to the Lagrange elements [33], and the spaces ■prA"(T) and 7'~A"(T) 
correspond to the space of discontinuous piecewise polynomials of degree < r and 
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< r — 1, respectively. When n = 2, the spaces 7'rA^(T) and V^h}{T) corre- 
spond to the Brezzi-Douglas-Marini iJ(div) elements of degree < r, introduced in 
[2F and the Raviart-Thomas 7? (div) elements of degree < r — 1 introduced in [75]. 
These spaces were generalized to three dimensions by Nedelec [75] , [73] . The spaces 
VrA^iT), VrA'^iT), P,7Ai(r), and V~K^{T) then correspond to the Nedelec 2nd 
kind _ff (curl) and iJ(div) elements of degree < r and the Nedelec 1st kind iJ (curl) 
and i7(div) elements of degree < r — 1, respectively. 

5.3. Computational bases. This subsection relates to the implementation of the 
V and families of finite element differential forms. It is not essential to the rest of 
the paper. Because the spaces V~A'^{T) and 7^,.A'^(T) were constructed through a 
finite element assembly procedure (shape functions and degrees of freedom) , we have 
at hand a basis for their dual spaces. Consider the space ■p^A'"'(T), for example, 
with reference to the decomposition of the dual space of the shape function space 



'P~A''{T) given in Theorem 5.5 Choose any / £ A(T) of dimension d > k, and 
choose any convenient basis for Vr+k-d-A'^^^iJ). For each of the basis functions ry, 
we obtain an element of V^K^{T)*: 

tr fLu AT], ujeV~K^{T). 



f 

(This is meaningful since the integral is, by construction, single- valued.) Taking 
the union over / e A(T) of the sets of elements of P~A''{T)* obtained in this 
way, gives a basis for that space. An interesting case is that of the Whitney forms 
Vi h^{T). In this case there is exactly one dual basis function for each / G Afc(T), 
namely w i— s- tr^ 

For computation with finite elements, we also need a basis for the finite element 
space itself, not only for the dual space. One possibility which is commonly used 
is to use the dual basis to the basis for the dual space just discussed. Given an 
element of the basis of the dual space, associated to some / € A(T), it is easy to 
check that the corresponding basis function of the finite element space vanishes on 
all simplices T £ T which do not contain /. Thus we have a local basis, which is 
very efficient for computation. 

In the case of the Whitney space 7^]~A'^(T), the dual basis can be written down 
very easily. We begin with the standard dual basis for ViA'^{Th): the piecewise 
linear function associated to the vertex Xi is determined by Xi{xj) = dij (after 
picking some numbering xi, . . . , xn of the vertices). Then given a fc-face / with 
vertices 2:^(0), • ■ • ^^^{k) we define the Whitney form [55] p- 229] 
fc 

" ' I ■ - rfA^(o) A • • • A dA<^(,) A • • • A dA^(fc). 



For higher degree finite element spaces it does not seem possible to write down 
the dual basis explicitly, and it must be computed. The computation comes down 
to inverting a matrix of size d x d where d is the dimension of the space of shape 
functions. This can be carried out once on a single reference simplex and the result 
transferred to any simplex via an affine transformation. 

An alternative to the dual basis, which is often preferred, is to use a basis for the 
finite element space which can be written explicitly in terms of barycentric coordi- 
nates. In particular, for the Lagrange finite element space ■PrA°(T) consisting of 
continuous piecewise polynomials of degree at most r, one often uses the Bernstein 
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basis, defined piecewise by monomials in the barycentric coordinates, instead of the 
dual or Lagrange basis. It turns out that explicit bases analogous to the Bernstein 
basis can be given for all the finite element spaces in the V and V~ families, as was 
shown in (lOj . Here we content ourselves with displaying a few typical cases. 
Bases for the spaces V^h}{T) and Vr^^{T) are summarized in Tables 



5.1 



and 

|5.2| respectively, for n = 3 dimensions and polynomial degrees r = 1, 2, and 3. To 
explain the presentation, we interpret the second line of Table [STT] We are assuming 
that T is a triangulation of a 3-dimensional polyhedron. The table indicates that 
for the space 7^^A^(T), there are two basis functions associated to each edge of 
T and two basis functions associated to each 2-dimensional face of T. If an edge 
has vertices Xi and Xj with i < j, then on an simplex T containing the edge, the 
corresponding basis functions are given by Xi4>ij and \j4>iji where and \j are 
the barycentric coordinates functions on T equal to 1 at the vertices Xi and Xj, 
respectively, and 

= d\j — Xj dXi 

is the Whitney form associated to the edge. Similarly, if T contains a face with 
vertices Xi, Xj, Xk, i < j < k, then on T the two basis functions associated to the 
face are given by \k4'ij and Xj4>ik- 

Table 5.1. Basis for the spaces V^A^, n = 3. 



r Edge [xi,Xj] Face [xi, xj , x^] Tet [xi, xj , Xf^, xi] 

1 4>ij 

3 {A?, A^, AiAj}(^ij {Ai, Aj , Afc}Afe0ij \i^\i4>ij, \j\i4>ik, 



{Xi, Xj , Xk}Xj4>ik 



^i^k4>il 



Table 5.2. Basis for the space VrA^, n = 3. 



Tet [xi,Xj,x^:,xi] 



Xi^dXiAdXj, XjdXiAdXi^, XidXj A dXf; 



X'^dXi A dXj, 
X^dXi A dAfe, 



XldX, 



XjXkdXi A (i(Afc - Xj) 
XiXjd{Xj - Xi) A dAfc 
AiAfedAj A d(Afc - A^) 

xU\ 



A dXi, 



A dXj, X^dXi A dX]^, ■\ia,Aj /\ aA^ 
X'^XkdX, A d{2Afe - Xj), XjXldXi A (i(Afc - 2Xj) 



XfXjd{2Xj 
X^X]d{Xj - 



-Ai)AdAfc, XfXkdXj A d{2Xk - Xi) 
2Xi)AdXk, XiXldXj A d(Xk ~ 2Xi) 
XiXjXkd(2Xj - A, - Afe) A d{2Xk - X^ - Xj) 



X^XidXiAdXj, XjXidXiAdX^. 
XjXj^dXi A dXi, XiXidXj A dXf^ 
XiXi^dXj A dXi, XiXjdXf; A dXi 

{Afe, AJAfcAjdAi A dXj 
{-^ji -^fci ^i}^j^idXi A dXk 
{-^ji ^kj ^i}^j^kii^i A dXi 
{Xi, Xj, Xk, Xi}XiXidXj A dXk 
{Ai, Aj , Afe, AJAiAfedAj A dXi 
{Xi,Xj,X^.,Xl}XiXjdX^: A dXi 
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5.4. Approximation properties. To apply the general theory for approximation 
of Hilbert complexes to the finite element exterior calculus, outlined above, we need 
to construct bounded cochain projections from L^A'^(ri) onto the various finite 



element spaces. The error estimates given in Theorems |3.9| and 3.11 bound the 
error in the finite element solution in terms of the approximation error afforded by 
the subspaces. In this subsection, we show that the spaces Vr^'^iT) and 'P~ K^{T) 
provide optimal order approximation of differential /c-forms as the mesh size tends 
to zero, where optimal order means that the rate of convergence obtained is the 
highest possible given the degree of the piecewise polynomials and the smoothness 
of the form being approximated. In the next subsection, we shall construct Lp'- 
bounded cochain projections which attain the same accuracy. 

Let {7^} be a family of simplicial triangulations of C M", indexed by decreasing 
values of the mesh parameter h given as h = maxj^gT-^ diamT. We will assume 
throughout that the family {7/i} is shape regular, i.e., that the ratio of the volume 
of the circumscribed to the inscribed ball associated to any element T is bounded 
uniformly for all the simplices in all the triangulations of the family. Parts of the 
construction below simplify if we assume, in addition, that the family {Th} is quasi- 
uniform, i.e., that the ratio /i/diamT is bounded for all T G T^, uniformly over 
the family. However, we do not require quasi-uniformity, only shape-regularity. 
Throughout this and the following subsection, h'^^ will denote a subspace of 

for each 



HK^ ^(fi) and a subspace of HA''{fl). Motivated by Section 5.1.4 
r = 0, 1, . . . , we consider the following possible pairs of spaces: 

(74) At' -I or f' 

[v-+,A'^-'{T)j [VrAHT) (ifr>0). 



Note that contains all polynomials of degree r, but not r-l- 1 in each case. As we 
have seen in Section |5.2[ the degrees of freedom define canonical projection opera- 
tors Ifi = iji mapping CA'^(ri) boundedly to A^ and commuting with the exterior 
derivative. However, the operators Ih do not extend boundedly to all of L^A'^(n), 
or even to HA''{il) and so do not fulfill the requirements of bounded cochain pro- 
jections of the abstract theory developed in Section [3] In this subsection, we will 
describe the Clement interpolant, which is a modification of the canonical projec- 
tion which is bounded on L^A'^ and can be used to establish the approximation 
properties of the finite element space. However, the Clement interpolant is not 
a projection and does not commute with the exterior derivative. So in the next 
subsection we construct a further modification which regains these properties and 
maintains the approximation properties of the Clement interpolant. 

The difficulty with the canonical interpolation operators Ih is that they make use 
of traces onto lower dimensional simplexes, and as a consequence, they cannot be 
extended boundedly to L^A''{il) (except for k = n). The Clement interpolant [M] 
is a classical tool of finite element theory developed to overcome this problem in the 
case of 0-forms. To define this operator (for general k), we need some additional 
notation. For any / e A{Th), we let fi/ C be the union of elements containing 



nf^\J{T\Te%,feA{T)}, 
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and let Pj : L'^K^{^f) Vr^^i^f) be the projection onto polynomial fc- forms 
of degree at most r. For uj G L^A''(ri), we determine the Clement interpolant 
IhOj :— ij^iv £ by specifying (j){Ihijj) for each degree of freedom of the space 



A^; see Section 5.2 Namely, if is a degree of freedom for A| associated to 
/ e A(7^), we take 

where ujf = i^\Q,f ■ 

The Clement interpolant is local in the sense that for any T £ T^, In'^W is 
determined by wj^*, where T* = 1J{^/ I / ^ A(^) }• fact, by scaling, as in [5J 
Section 5.3], for example, it can be seen that 

ll^/i'^llL2Afc(T) < Co ^ ||Pht^||L2Afc(n^) < Ci||cj||i2Afc(T*)j 
/eA(T) 

where the constants co,ci may depend on the polynomial degree r and the di- 
mension n, but are independent of T G 7/j and h, thanks to the shape regu- 
larity assumption. Therefore, the Clement interpolant is uniformly bounded in 

ll^/iW||L2Afc(o) < C2||w||L2Afc(o) 

(with the constant independent of h). 

Another key property of the Clement interpolant is that it preserves polynomials 
locally in the sense that //jw|t = w|t if w € VrA^iT*). It is then standard, using 
the Bramble-Hilbert lemma [23] and scaling, to obtain the following error estimate. 

Theorem 5.8. Assume that Af^ is either V^^^A^^ {%) , or, if r > 1, VrA^iTh). 
Then there is a constant c, independent of h, such that the Clement interpolant 
ijl : L'^A^{Q) A^ satisfies the bound 

||cj - /^cj||i2A'.(n) < c/i"|cj|H=A'=(n), ^ e H'A''{n), 
for < s <r + 1. 

Note that the estimate implies that any sufficiently smooth /c-form is approx- 
imated by elements A^ with order 0(ft,''+^) in L^. Since the polynomial spaces 
used to construct A^^ contain VrA'' but not Vr+iA'', this is the optimal order of 
approximation . 

5.5. Bounded cochain projections. The Clement interpolant Ih is both uni- 
formly bounded in and gives optimal error bounds for smooth functions. How- 
ever, it is not a bounded cochain projection in the sense of the theory of Section 
3. Indeed, it is neither a projection operator — it does not leave A^ invariant — nor 
does it commute with the exterior derivative. Therefore, to construct bounded 
cochain projections, we consider another modification of the canonical projection 
III, in which the operator 1^ is combined with a smoothing operator. 

This construction — key ingredients of which were contributed by Schoberl [73] 
and Christiansen |31j — was discussed in detail in [5^ (where it was called a smoothed 
projection), under the additional assumption that the family of triangulations {Th} 
is quasi-uniform, and then in [32] in the general shape regular case. Therefore, we 
will just give a brief outline of this construction here. 

Let A^ be one of the spaces 'P~^-i^A''{Th) or, if r > 1, VrA^{Th), and let Ih ■ 
CA'"'(ri) A{^ be the corresponding canonical interpolant. To define an appropriate 
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smoothing operator, we let p : M" — )■ M be a nonnegative smooth function with 
support in the unit ball and with integral equal to one. In the quasi-uniform case, 
we utilize a standard convolution operator of the form 

I p{y)^{x + 5y) dy, 
"'|y|<i 

mapping L?'A^{Cl) into C°°h^{^), where the domain 17 D is sufficiently large 
so that the convolution operator is well defined. We set the smoothing parameter 
5 = eh, where the proportionality constant e > is a parameter to be chosen. In 
the more general shape regular case, we need to generalize this operator slightly. 
As in [32] . we introduce a Lipschitz continuous function : — > M+, with Lips- 
chitz constant bounded uniformly in h, such that gh{x) is uniformly equivalent to 
diamT for T d and a; £ T. Hence, ghlr approximates diamT. The appropriate 
smoothing operator in the general case, mapping L^A'^(fi) into CA'''(ri), is now 
given by 

uj^ f piym^yycohdy, 

where the map : E" M" is defined by 

^ x + egh{x)y. 

Since this smoothing operator is defined as an average of puUbacks, it will indeed 
commute with the exterior derivative. By combining it with an appropriate exten- 
sion operator E : Lp'K^{Q) L^A'^(f2), constructed such that it commutes with the 
exterior derivative, we obtain an operator i?^ : L^A'^(ri) CA'^(ri). The operator 
= Ih o Rj^ maps L^A'^{il) into A^, and commutes with the exterior derivative. 
Furthermore, for each e > 0, the operators are bounded in L^, uniformly in h. 

However, the operator is not invariant on the subspace A^. The remedy to 
fix this is to establish that the operators Q^|a^ ■ -^h ^ -^h converge to the identity 
in the operator norm as e tends to zero, uniformly in h. Therefore, for a fixed 
e > 0, taken sufficiently small, this operator has an inverse : A^ ^ A^. The 
desired projection operator, tt/j — n^, is now given ds tt^ — Jf^o Q^. The operators 
TT/i are uniformly bounded with respect to h as operators in £(L^A'^(ri), L?K^{^)). 
Furthermore, since they are projections onto A^, we obtain 

(75) ||w - 7rftw||i2Afc(n) < inf ||(/-7r/i)(w-/i)||i2Afc(j2). 

Based on these considerations, we obtain the the following theorem. Cf. ^ Theorem 
5.6] and [32l Corollary 5.3]. 

Theorem 5.9. 1. Let A^ he one of the spaces V^^+^A''' {%) or, ifr>\, VrA'^iTh) 
and TT^ : _L^A'^'(il) — > A^ the smoothed projection operator constructed above. Then 
TT^ is a projection onto A^ and satisfies 

- 4^^\\L^^''{n} < ch'\\uj\\HsAk(n), ^ e H'K^{n), 

for < s < r + 1. Moreover, for all lu G L^A'^(ri), Tr^tj —> uj in L'^ as h ^ 0. 
2. Let A^ he one of the spaces VrA''{Th) or 'P;:A''{Th) with r > 1. Then 

\\d{uj - 7r,';w)||i2A'o(o) < ch'\\duj\\H^AHn). ^ e H'A''{n), 
for < s < r. 
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3. Let ^ and be taken as in (74 1 and tt^ ^ and tt^ the corresponding 



smoothed projections. Then dn. 



k-l 



Proof. We have already established that the tt^ are uniformly bounded projections 
which commute with d. The error estimate in the first statement of the theorem 
then follows from (751. Statement 2 follows from 1 and 3. 



□ 



5.6. Approximation of the de Rham complex and the Hodge Laplacian. 

Let AjJ be the space of piecewise polynomial n-forms of degree at most s for some 



nonnegative integer s. 
de Rham complexes 



We showed in Section 5.2 that there are 2" distinct discrete 



(76) 







AO 



A" 



A" 



0, 



with each of the mapping Aj'j^^ — > A^ being of the form (73 1 for an appropriate 
choice of s. Making use of the bounded projections, we obtain in each case a 
commuting diagram 







K 



iJA"(r2) ^ 



0. 



Thus the projections give a cochain projection from the de Rham complex to 
the discrete de Rham complex, and so induce a surjection on cohomology. Since 
lim/,^0 \\(^-T^hw\\HA''{n) 



3.4 



(with V = HA and Vh = Ah) 



0, we can use Theorem i 
to see immediately that for h sufficiently small, this is in each case an isomorphism 
on cohomology. 

The simplest finite element de Rham complex is the complex of Whitney forms. 



(77) 



O^V^A°iTh)^r^A\%) 



V^A^\%) ^ 0. 



Via the basic Whitney forms cj)f, the space Vi A^{Th) is isomorphic to the space 
of simplicial cochains of dimension k associated to the triangulation 7/^. In this 



way, the complex (77 1 of Whitney forms is isomorphic, as a cochain complex, to 



the simplicial cochain complex, and its cohomology is isomorphic to the simplicial 
cohomology of Q.. Thus Theorem |3.4| when applied to the Whitney forms, gives 
an isomorphism between the de Rham cohomology and simplicial cohomology for 
a sufficiently fine triangulation. Since the simplicial cohomology is independent of 
the triangulation (being isomorphic to the singular cohomology), for the complex 
of Whitney forms, the isomorphism on cohomology given by Theorem |3 .4| holds for 
any triangulation, not just one sufficiently fine. This proves de Rham's theorem on 
the equality of de Rham and simplicial cohomology. 

For all of the discrete complexes ([76|, as for the Whitney forms, the cohomology 
is independent of the triangulation (and equal to the de Rham cohomology). To see 



this, following |31j . we consider the Whitney forms complex (77) as a subcomplex 



of (76). The canonical projections Ih define cochain projections. Note that the Ih 
are defined on the finite element spaces A^, because all of the trace moments they 
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require are single- valued on Ajj. From the commuting diagram 



ih 



ih 



Q-.V^A\Th) V^A\Th) ••• — ^ PrA"(r;,)->0, 

we conclude that the cohomology of the top row, which we have already seen to be 
an image of the de Rham cohomology, maps onto the cohomology of the bottom 
row, which is isomorphic to the de Rham cohomology. Hence the dimension of 
all the corresponding cohomology groups are equal and both cochain projections 
induce an isomorphism on cohomology. 



Consider now the numerical solution of the Hodge Laplacian problem ( 59 ) . We 
approximate this by a Galerkin method: Find (Jh G ^/J^^i Uh G A^^, and ph £ S)^ 
such that 

{cjh,T) - {dT.Uh) reA{;"\ 
{dah,v) + {duh,dv) + {v,ph) = (/,«), v G AJ;, 



h- 



For the finite element spaces, we may choose A^ ^ and Ajj to be any of the pairs 
given in (74). We have verified that d maps A^^^ into A^ and this map can be 



extended to a full subcomplex of the de Rham complex admitting a bounded cochain 
projection. Therefore we may combine the error estimates derived in the abstract 
setting in Sections |3.5| and |3.6| and approximation estimates from Section |5.5| to 
obtain convergence and rates of convergence. The rates of convergence are limited 
by three factors: (1) the smoothness of the data /, (2) the amount of elliptic 
regularity (determined by the smoothness of the domain), and (3) the degree of 
complete polynomials contained in the finite element shape functions. Specifically, 
from Theorem |3 . 1 1| we get the following result. Assume the regularity estimate 

||M|j_H-t+2(o) + \\p\\Ht+^(n) + + lk||//t+i(o) + \\do\\Hi(n) < C'||/||_ff*(n), 

holds for < < < tmax- Then for < s < i^ax the following estimates hold: 

\\d{a-ah)\\ <C/in|/||H.(o), if.s<r + l, 




s<r, Al^V;^,A^{T), 
s<r-l, Al^VrA^{T), 



\\u-Uh\\ + \\p-Ph\\ < 



ChWfW, Al^V^A^T), 
C/i'*+^||/||//s(Q) if s < r — 1, otherwise. 

Thus, the error in each case is the optimal order allowed by the subspace if we have 
sufficient elliptic regularity. 

Concerning the eigenvalue problem, applying Theorem |3.21| we immediately 
obtain the following error estimate: 
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This estimate does not apply in case — A''{T), i.e., the Whitney forms. In 
that case, we get the estimate: 



|A-A„| < 



Ch\\uf, 

Ch^u\\\\u\\Hiin)- 



6. Variations of the de Rham complex 

In Sections |4] and |5] we apphed the abstract theory of Section [3] to the Hilbert 
complex obtained by choosing W'' to be L^A^-^^) ^^^d d*"' to be the exterior de- 
rivative with domain V'^ = HA^{^1). This led to finite element approximations of 
certain boundary value problems for the Hodge Laplacian, and a variety of related 
problems (03 and S* problems, eigenvalue problems). In this section, we show 
how some small variations of these choices lead to methods for important related 
problems, namely generalizations of the Hodge Laplacian with variable coefficients 
and problems with essential boundary conditions. 

6.1. Variable coefRcients. In this section, d'^ is again taken to be the exterior 
derivative with domain V'^ — HA'^{n) and W'' is again taken to be the space L^h^ , 
but is furnished with an inner product which is equivalent to, but not equal to, 
the standard L^A*^ inner product. Specifically, we let : be a bounded, 

symmetric, positive definite operator and choose 

{u,v)^/k = {a^u,v)L2fi^k. 

With this choice, we can check that d^. — {d''~^)* is given by 

dlu={a''-Y'Ska''u, 

with the domain of d^. given by 

= {ue L^A'^ifl) I a'^u £ H*A''{n) } = {a'')~^H*A''{n). 

The harmonic forms are then 

^l''^ - S'' n 3fc - { w e HA'' n {a'')-^H*A'' \du^O, dal^u = }. 

The Hodge Laplacian d*du + dd*u is now realized as 

d{a''~^)-^5a''u + {a'')-^5a''+^du. 

As an example of the utility of this generalization, we consider the S| eigenvalue 
problem in 3 dimensions. We take — e and a? = where e and ^ are symmetric 
positive definite 3x3 matrix fields. Then the 23^; eigenvalue problem is 

curl/i^^ curlii = \eu. 

This is the standard Maxwell eigenvalue problem where e is the dielectric tensor of 
the material and /i is its magnetic permeability. These are scalar for an isotropic 
material, but given by matrices in general. They are constant for a homogeneous 
material, but vary if the material is not constant across the domain. 

We notice that the Hilbert complex in this case is exactly the de Rham complex 



(52 1, except that the spaces are furnished with different, but equivalent, inner prod- 
ucts. Thus the finite element de Rham subcomplexes we constructed in Section [5] 
apply for this problem, and satisfy the subcomplex and bounded cochain projection 
properties. In other words, the same finite element spaces which were developed for 
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the constant coefficient Hodge Laplacian, can be applied equally well for variable 
coefficient problems. 

6.2. The de Rham complex with boundary conditions. For another impor- 
tant application of the abstract theory of Section|3] we again choose = L^h^in) 
and to be the exterior derivative, but this time we take its domain to be 
the space HA''{il). This space is dense in W'^, and is a closed subspace of HA'', 
so is complete, and hence d is again a closed operator with this domain. Thus we 
obtain the following complex of Hilbert spaces, called the de Rham complex with 
boundary conditions: 

(78) 0-^ HA°{n) ^ HA\n) ^ ■■■ ^ HA"{n) 



0. 



In Section 



4.2 



we used the fact that HA'' n H*A'' is compactly included in L^A''. 
Replacing fc by n — fc and applying the Hodge star operator, we conclude the 
analogous result, that HA'' n H*A'' is compactly included in L^A''. Then, just 
as in Section 4.2 all the results of Section |3] apply to the de Rham complex with 
boundary conditions. In particular, we have the Hodge decomposition of L^A'^ and 
of HA'', the Poincare inequality, well-posedness of the mixed formulation of the 



Hodge Laplacian, and all the approximation results established in Sections 3.3 3.6 



We denote the spaces of fc-cocycles, fc-coboundaries, and harmonic fc- forms as 
^'' = {uje HA'' \duj^O}, S'^' = { dr; I ?7 e HA''-^ }, 

f)''^{uje5''\{oj,fi) = 0, ^ie%''}. 



Now, from (53 1, we have that 
(79) 
Therefore 

%''^l2 --{uje L^A'' I {uj,d?]) 



Vr/ e HA''^^ } = { e H*A'' \5oj = Q}. 



This gives a concrete characterization of the harmonic forms analogous to ( 58 1 



(80) 



{uj e HA''{n) n H*A''{n) \duj = o,SLo^ o}. 



Now, from (54 1 
with duj 



( [55| , and (50 1, the Hodge star operator maps the set of w G i?A"~'^ 
isomorphic ally onto the set of /z G H*A'' with Sjj, — 0, and maps the 
set of G 7i*A"^'^ with 6u! — isomorphically onto the set of /i G HA'' with 
dfi ~ 0. Comparing (58 1 and (80 1, we conclude that the Hodge star operator maps 
Sj"^'' isomorphically onto S)'' 



particular, dim^'' is the (n 



This isomorphism is called Poincare duality. 
fc)th Betti number of il. 



In 



Finally, from ( 79 1 , we can characterize d* , the adjoint of the exterior derivative 
d viewed as an unbounded operator L^A''^^ L^A'' with domain HA''. Namely, 



in this case d* has domain H*A'' and coincides with the operator S defined in (50 1. 
By contrast, when we took the domain of d to be all of HA'', the domain of d* 
turned out to be the smaller space H* A'' . 

We now consider the Hodge Laplacian problem given in the abstract case by 
(18 1. In this case, we get that {a,u,p) G HA''^^ x HA'' x ^'^ is a solution if and 



only if 

(81) 

(82) 



: 5u, da - 
tva = 0, 



- Sdu = f 
tr w = 



- p in Q, 
on dH., 
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(83) 



Note that now both boundary conditions are essential in the sense that they are 
implied by membership in the spaces HA^~^ and HA^ where the solution is sought, 
rather than by the variational formulation. 

To make things more concrete, we restrict to a domain C M'^, and consider 
the Hodge Laplacian for fc- forms, /c 0, 1,2, and 3. We also discuss the 5B* and S 



problems given by (19 1 and (20 1 for each k. Again we have 



(f' ~ grad, d} ~ curl, d? 
and the de Rham complex ([78| is now realized as 



div, 



H^{n) 



H{div;n) ^ L'^{n) -> 0, 



where 

H(curl; ft) = {u : ft ^ M.^ \ u e L'^ , cmiu € u x n = on dn }, 
H{div;n) ^ {u -.n -^R^lu e L"^, divueL^, u ■ n ^ on dn}. 
6.2.1. The Hodge Laplacian for k = 0. For fc = 0, the boundary value problem 



(81l-(83) is the Dirichlet problem for the ordinary scalar Laplacian. The space 



HA ^ is understood to be 0, so a vanishes. The space Sj^ of harmonic 0-forms 



vanishes. The first differential equation of (81 1 vanishes, and the second gives 
Poisson's equation 

— div grad u = f in f2. 



Similarly, the first boundary condition in (82 1 vanishes, while the second is the 
Dirichlet condition li on dil. The side condition (83 1 is automatically satis- 



fied. Nothing additional is obtained by considering the split into the 5B* and ?B 
subproblems, since the latter is trivial. 



(84) 



6.2.2. The Hodge Laplacian for k — I. In this case the differential equations and 
boundary conditions are 

(7 — ~ div u, grad a + curl curl u ~ f — p in fi, 
(7 = 0, wxn = on 917, 

which is a formulation of a boundary value problem for the vector Laplacian 
curl curl — grad div. The solution is determined uniquely by the additional con- 
dition that u be orthogonal to f)^ = io^, i.e., those vector fields satisfying 

curlp = 0, divp = in f2, p x n = 0, on d^l. 



The 58^ problem (19) is defined for L^ vector fields / which are orthogonal to 
both gradients of functions in H^{il) and the vector fields in Sj^. In that case, the 
solution to (84) has <t — and p = 0, while u satisfies 

curl curl u = /, divM = in fi, u x n = on dfl. 

The orthogonality condition u J- Sj^ again determines the solution uniquely. 

Next we turn to the problem. For source functions of the form / — grad F for 



some F g 77 ^ , ( 84 1 reduces to the problem of finding a £ H^ and u G = grad H^ 
such that: 



a — — div u, grad cr = f in fi. 
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The differential equations may be simplified to — grad div u ~ f and the condition 
that u e 58^ can be replaced by the differential equation curlu = and the bound- 
ary condition u x n = 0, together with orthogonality to Sj^. Now grad(cr — F) = 
and a — F = on dfl, so a = F. Hence, we may rewrite the system as 

— div u — F, curlu = in f2, u x n = 0, on dfl, 

which, again, has a unique solution subject to orthogonality to 

6.2.3. The Hodge Laplacian for k ~ 2. The differential equations and boundary 
conditions are 

a — curl u, curl a — grad div u — f ^ p in fi, 
axn = 0, u-n = on dil. 
This is again a formulation of a boundary value problem for the vector Laplacian 



curl curl — grad div, but with different boundary conditions than for (84 1, and this 
time stated in terms of two vector variables, rather than one vector and one scalar. 
This time uniqueness is obtained by imposing orthogonality to = i^^, the space 
of vector fields satisfying 

curlp — 0, div p =0 in fi, p ■ n = 0, on dfl. 

The problem arises for source functions of the form / = grad F for some 
F G which we may take to have integral 0. We find a — 0, and u solves 

— div u — F, curl u — in fl, u • n = 0, on dfl, 

i.e., the same differential equation as for *8^, but with different boundary conditions 
and, of course, now uniqueness is determined by orthogonality to Sj^. 

ff div / = 0, / • n = on dH,, and / _L S^^, we get the 05^ problem for which the 
differential equations are a = curlu, curler — f and the condition divu = arising 
from the membership of u in 05^. Thus u solves 

curl curl u ~ /, divu = in 51, it • n = 0, on 9f2, 

the same differential equation as for 03^, but with different boundary conditions. 

6.2.4. The Hodge Laplacian for k ^ 3. The space S)'^ consists of the constants (or, if 
il is not connected, of functions which are constant on each connected component). 
The Hodge Laplacian problem is 

cr = — grad u, div cr = / — p in il, a-n — Oondfl, / uqdx^O, q £ Sj^ , 

Jn 

which is the Neumann problem for Poisson's equation, where the (piecewise) con- 
stant p is required for there to exist a solution and the final condition fixes a unique 
solution. 

7. The elasticity complex 

fn this section, we present another complex, which we call the elasticity complex, 
and apply it to the development of numerical methods for linear elasticity, fn 
contrast to the other examples presented above, the elasticity complex is not a 
simple variant of the de Rham complex, and, in particular, it involves a differential 
operator of second order. However, the two complexes are related in a subtle 
manner, via a construction known as the Bernstein-Bernstein-Gelfand resolution, 
explained for example in [5J IHl 113] and references given there. 
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The equations of elasticity are of great importance in modeling solid structures, 
and the need to solve them in engineering applications was probably the primary 
reason for the development of the finite element method in the 1960s. The sim- 
plest finite element methods for elasticity are based on displacement approaches, in 
which the displacement vector field is characterized as the minimizer of an elastic 
energy functional. The design and analysis of displacement finite element meth- 
ods is standard and discussed in many textbooks, (e.g., [33i). However, for more 
general material models, arising, for example, for incompressible materials or some 
viscoelastic or plastic materials, the displacement method is either not feasible or 
performs poorly, and a mixed approach, in which the elastic stress and the dis- 
placement are taken as unknowns, is the natural alternative. In fact, mixed finite 
element methods for elasticity were proposed in the earliest days of finite elements 
|51j . and stable finite elements for the mixed formulation of elasticity have been 
sought for over four decades. These proved very elusive. Indeed, one of the motiva- 
tions of the pioneering work of Raviart and Thomas '76j on mixed finite elements 
for the Laplacian, was the hope that the solution to this easier problem would pave 
the way to such elements for elasticity, and there were many attempts to generalize 
their elements to the elasticity system [4J [SJ [HOI [HI IHll- However, it was not 
until 2002 that the first stable mixed finite elements for elasticity using polynomial 
shape functions were discovered by two of the present authors |TT], based on tech- 
niques and insights from the then nascent finite element exterior calculus. See also 
[7]. They developed and analyzed a family of methods for plane (2-dimensional) 
elasticity. The basic elements in the family involve a space of shape functions inter- 
mediate between quadratics and cubics for the stress, with 24 degrees of freedom, 
while the displacement field consists of the complete space of piecewise linear vector 
fields, without continuity constraints. It is these elements we will describe here. 

The plane elasticity elements of |11) were generalized to three dimensions in 
[3]. See also [1]. These elements are, however, quite complicated (162 degrees of 
freedom per tetrahedron for the stress) . A much more promising approach involves a 
variant Hilbert complex for elasticity, the elasticity complex with weak symmetry. 
In [5] and [S], stable families of finite elements are derived for this formulation 
in n-dimensions, again using techniques from the finite element exterior calculus. 
For the lowest order stable elements, the stress space is discretized by piecewise 
linears with just 12 degrees of freedom per triangle in two dimensions, 36 per 
tetrahedron in three dimensions, while the other variables are approximated by 
piecewise constants. 

7.1. The elasticity system. The equations of linear elasticity arise as a system 
consisting of a constitutive equation and an equilibrium equation: 

(85) Aa — eu, diva^/ in 17. 

Here the unknowns a and u denote the stress and displacement fields engendered by 
a body force / acting on a linearly elastic body which occupies a region C K". The 
constitutive equation posits a linear relationship between the linearized deformation 
or strains due to the displacement and the stresses. The equilibrium equation states 
that the stresses, which measure the internal forces in the body, are in equilibrium 
with the externally applied force. The stress field a takes values in the space S := 
M"y^ of symmetric matrices and the displacement field u takes values in V := K". 
The differential operator e is the symmetric part of the gradient, the div operator 
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is applied row-wise to a matrix, and the compliance tensor A — A(x) : § — )■ S is a 
bounded and symmetric, uniformly positive definite operator reflecting the material 
properties of the body. If the body is clamped on the boundary dfl of fi, then the 



appropriate boundary condition for the system (85 1 is u = on dfl. For simplicity, 
this boundary condition will be assumed here. 

The pair (cr, u) can alternatively be characterized as the unique critical point of 
the Hellinger-Reissner functional 

/(t, v) — / {-p:AT : T + div t ■ v ~ f ■ v) dx. 

The critical point is sought among r S H{div; il; S), the space of square-integrable 
symmetric matrix fields with square-integrable divergence, and v G L^(f2;V), the 
space of square-integrable vector fields. Equivalently, the pair (cr, u) € H{div;^l;S) 



xL^(ri;V) is the unique solution to the weak formulation of the system (85 I: 
J^{A(j : T + div T ■ u) dx = 0, r e iJ(div; fi; §), 

(86) 

J^^diva-vdx = J^^f-vdx, v L'^{fl;W). 

7.2. The elasticity complex and its discretization. The discretization of the 



problem (86 1 is closely tied to discretization of the elasticity complex. For a two 
dimensional domain fi, the elasticity complex takes the form 

(87) O^H^{n) — ^ i/(div;r!;§) L^{n;R^)^0. 

The operator J is a second-order differential operator mapping scalar fields into 
symmetric matrix fields, namely the rotated Hessian 

""{-Sk, -|'")=ofe'«is»d,)o- 0.(_°j J). 

We now show that this complex is closed, i.e., that the range of the two operators 
J and div are closed subspaces. The null space of J is just the space Vi of linear 
polynomials and we have ||g||//2 < c|| J(7||i2 for all q S H^{U) with q J-Vi, which 
implies that the range of J is closed. We now show that div maps J?(div;ri;§) 
onto L^(n;]R^), using the fact that the scalar- valued divergence maps iJ^(ri;IR^) 
onto L^(ri). Given / G L'^(ri;M^) the latter result implies the existence of r S 
i?"'^(f2; M^^^) with divr = /. This t may not be symmetric. However, for any 
V e -ff^(ri;]R^), we let curlt; denote the matrix with ith row equal to curlw; = 
{—dvi/dx2,dvi/dxi). Then div(r -I- curlti) = / and r-|-curlw € H{div;U;E>) if 
divw = T21 — Ti2. Thus div : H{d\v; L^(ri;K^) is indeed surjective. Setting 

W° = L^n), = L^{^;S), and W"^ = L'^{n;R^), we obtain a closed Hilbert 



complex with domain complex ( 87 1 



Now we describe only the simplest discretization from llj. The space V^^ C 
F° = iJ2(f]) is the finite element space whose shape functions on each triangle T 
are the quintic polynomials 'Pz{T) and whose degrees of freedom are the values, 
first, and second derivatives at each vertex, together with the mean value of the 
normal derivatives along each edge, which gives 21 = dim7-'5(T) degrees of freedom 
per triangle. The resulting finite element space is referred to in the finite element 
literature as the Argyris space, and is easily seen to be contained in C^{Q), and 
hence in H^{yi). In fact, the Argyris space is in a certain sense the simplest finite 
element space contained in C^{yi). However the Argyris space does not coincide 
with the space of all piecewise quintics (which is not a finite element space, as 
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it cannot be defined via degrees of freedom) . Elements of the Argyris space have 
extra smoothness at the vertices, where they are C^. 

The space C = H{div;n;E>) consists of those continuous piecewise cu- 
bic symmetric matrix fields whose divergence is piecewise linear. The degrees of 
freedom determining an element a e Vfl are 

• the value of cr at each vertex, 

• the moments of degree and 1 of the two normal components of a on each 
edge, and 

• the moments of degree of cr on each triangle. 

Hence, the restriction of cr to a triangle T e 7^ is uniquely determined by 24 degrees 
of freedom. Finally, the space Vj^ C V'^ — L^(ri;M^) consists of all piecewise linear 
vector fields with respect to T^, without any imposed continuity. For degrees of 
freedom, we choose the moments of degrees and 1 for each component on each 
triangle T. Using the degrees of freedom, we also define projections ij^v G 
defined for sufliciently smooth v £ by (f>{ljlv) — (f){v) for all degrees of freedom 



It is straightforward to check that we have the subcomplex property jy,° C V/l 

ijl give a commuting 



and div V/l C and that the interpolation operators 1^ 



diagram 



nc" 



div 



div 



However, since l)^ involves point values of the second derivative and 1\ involves 
point values, these operators do not extend to bounded operators on the Hilbert 
spaces and . 

In order to apply the general theory of approximation of Hilbert complexes 
developed in Section [sj we will modify the interpolation operators in ([88 1 to 
obtain bounded cochain projections tt^ : for A: = 0, 1, 2. As in Section 5.5 

the TT^ will be obtained from via smoothing. 

7.3. Bounded cochain projections for the elasticity complex. First, we de- 
fine an appropriate puUback of an affine map, and show that it defines a cochain map 
for the elasticity sequence. Let : ^ fi' C be an affine map, so that DF = 
DF{x) is a constant 2x2 matrix. Let B be the matrix B = det{DF){DF)^^ = 
0{DFyO^ and define puUbacks F* = F* : W''{n') VF*=(f]) by 

F*v^voF, v£W°{fl'), 

F*v^ B{voF)B'^ , v£W^{n'), 

F*v = det{DF)B{vo F), veW^{n'). 

It is straightforward to check that, as long as F is affine, F* commutes with the 
differential operators J and div. In other words, we obtain the commuting diagram 



0-^^0(r!') 



v^{n') ^ 



"1 
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To define a smoothing operator which maps W'' into itself, we need to make 
the assumption that the domain is star-shaped with respect to some point in its 
interior (which we assume, without further loss of generality, to be the origin). 
We then let a = m&Xx^dn and define the dilation : by Fs{x) — 

x/{l + aS), S > 0. Then Fg maps the (5-neighborhood of fi, ils, into and so 
maps L^(ri) to L'^{ils). Composing with a standard mollification 

p{z)v{x + 6z) dz, 

\z\<l 

we obtain a smoothing operator Rs mapping W'^ into W'^. Here p : ^ E is 
a smooth, nonnegative function supported in the unit ball and with integral equal 
to one. This smoothing commutes with any differential operator with constant 
coefficients, in particular with J and div. 

The construction now proceeds as in Section [5.5| and we just give an outline of 
it. For simplicity, we assume a quasi-uniform family of meshes, although extension 



to general shape regular meshes could be made as in Section 5.5 Let e > be fixed. 
Define operators : W'' by Q'l = o R^^. The operator Q\ maps 

into . In fact, it can be seen by a scaling argument, that the operators Q\ are 
bounded as operators in C{W^,W^), uniformly in h. Furthermore, commutes 
with J and div, i.e.. 



(89) Qh° J ^ J °Qh and o div div o Qj^. 

However, the operators are not projections, since they do not restrict to the 
identity on V^. Therefore, we define the desired projections tt^ : W'^ — > V/^ as 
TT^ = (Q||y^fc)~^(5{;. To justify this definition, we need to argue that {Qh\v,'')~^ 
exists. In fact, for e sufficiently small, but fixed, {Qf^\yk)~^ is even uniformly 
bounded in C(Wl^,Wl^), where Wj^ — V/^, but equipped with the norm of W''. 
This follows from the fact that the operators Qf^\v^ converge to the identity with 
e in C{W^, Wjl), uniformly in h, estabhshed by a scaling argument. 

We conclude that if the parameter e > is fixed, but sufficiently small, then 
the operators tt^ are well-defined and uniformly bounded with respect to h in 
C{W^,W^). Furthermore, the diagram 

"4 "4 "4 



commutes, as a consequence of (89 1 



Thus, under the assumption of a star-shaped domain, we have established the 
existence of bounded cochain projections and so the results of Section [3^ apply to 
the elasticity complex. This means that the corresponding Hodge Laplacians are 
well-posed and that the discretizations are stable. In particular, using the spaces 
Vf^ and V^, we obtain a stable, convergent mixed discretization of the elasticity 
equations. 

Remark. Stability and convergence of the mixed finite element method for elasticity 
using the spaces Vf^ and described above is proven in [TT] without restricting to 
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star-shaped domains. However we do not know how to construct bounded cochain 
projections without this restriction. 
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